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Abstract 

A  high-order  discontinuous  Galerkin  finite  element  discretization  and  p-multigrid  solution 
procedure  for  the  compressible  Navier-Stokes  equations  are  presented.  The  discretization 
has  an  element-compact  stencil  such  that  only  elements  sharing  a  face  are  coupled,  regardless 
of  the  solution  space.  This  limited  coupling  maximizes  the  effectiveness  of  the  p- multigrid 
solver,  which  relies  on  an  element-line  Jacobi  smoother.  The  element-line  Jacobi  smoother 
solves  implicitly  on  lines  of  elements  formed  based  on  the  coupling  between  elements  in 
a  p  =  0  discretization  of  the  scalar  transport  equation.  Fourier  analysis  of  2-D  scalar 
convection-diffusion  shows  that  the  element-line  Jacobi  smoother  as  well  as  the  simpler 
element  Jacobi  smoother  are  stable  independent  of  p  and  flow  condition.  Mesh  refinement 
studies  for  simple  problems  with  analytic  solutions  demonstrate  that  the  discretization 
achieves  optimal  order  of  accuracy  of  0(hp+1).  A  subsonic,  airfoil  test  case  shows  that 
the  multigrid  convergence  rate  is  independent  of  p  but  weakly  dependent  on  h.  Finally, 
higher-order  is  shown  to  outperform  grid  refinement  in  terms  of  the  time  required  to  reach 
a  desired  accuracy  level. 
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Chapter  1 


Introduction 


The  simulation  of  complex  physical  phenomena  using  numerical  methods  has  become  an  in¬ 
valuable  part  of  modern  science  and  engineering.  The  utility  of  these  simulations  has  led  to 
the  evolution  of  ever  more  efficient  and  accurate  methods.  In  particular,  numerous  research 
efforts  have  been  aimed  at  developing  high-order  accurate  algorithms  for  solving  partial 
differential  equations.  These  efforts  have  led  to  many  types  of  numerical  schemes,  including 
higher-order  finite  difference  [33,  59,  56],  finite  volume  [7,  57],  and  finite  element  [8,  20,  6,  40] 
methods  for  both  structured  and  unstructured  meshes.  Despite  these  developments,  in  ap¬ 
plied  aerodynamics,  most  computational  fluid  dynamics  (CFD)  calculations  are  performed 
using  methods  that  are  at  best  second-order  accurate.  These  methods  are  very  costly,  both 
in  terms  of  computational  resources  and  time  required  to  reach  engineering-required  ac¬ 
curacy.  Higher-order  methods  are  of  interest  because  they  have  the  potential  to  provide 
significant  reductions  in  the  time  necessary  to  obtain  accurate  solutions.  Motivated  by 
this  potential,  the  goal  of  this  thesis  is  to  contribute  to  the  development  of  a  higher-order 
CFD  algorithm  which  is  practical  for  use  in  an  applied  aerodynamics  setting.  Specifically, 
the  thesis  details  a  high-order  discontinuous  Galerkin  (DG)  discretization  of  the  compress¬ 
ible  Navier-Stokes  equations  and  a  multigrid  solution  procedure  for  the  resulting  nonlinear 
discrete  system. 

1.1  Motivation 

While  CFD  has  matured  significantly  in  past  decades,  in  terms  of  time  and  computational 
resources,  large  aerodynamic  simulations  of  aerospace  vehicles  are  still  very  expensive.  In 
this  applied  aerodynamics  context,  the  discretization  of  the  Euler  and/or  Navier-Stokes 
equations  is  performed  almost  exclusively  by  finite  volume  methods.  The  evolution  of  these 
methods,  including  the  incorporation  of  upwinding  mechanisms  [51,  46,  52,  47,  53]  and 
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advances  in  solution  techniques  for  viscous  flows  [4,  41,  37,  38],  has  made  the  simulation  of 
complex  problems  possible.  However,  the  standard  algorithms  remain  at  best  second-order 
accurate,  meaning  that  the  error  decreases  as  0(h2). 

Moreover,  while  these  methods  are  used  heavily  in  aerospace  design  today,  the  time 
required  to  obtain  reliably  accurate  solutions  has  hindered  the  realization  of  the  full  potential 
of  CFD  in  the  design  process.  In  fact,  it  is  unclear  if  the  accuracy  of  current  second-order 
finite  volume  methods  is  sufficient  for  engineering  purposes.  The  results  of  the  two  AIAA 
Drag  Prediction  Workshops  (DPW)  [35,  31]  suggest  that  the  CFD  technology  in  use  today 
may  not  produce  adequate  accuracy  given  current  grids.  Numerous  authors  [54,  32,  22]  have 
shown  that  the  spread  of  the  drag  results  obtained  by  the  DPW  participants  is  unacceptable 
given  the  stringent  accuracy  requirements  of  aircraft  design. 

This  problem  could  be  alleviated  by  the  development  of  a  high-order  CFD  algorithm. 
Specifically,  a  high-order  method  could  reduce  the  gridding  requirements  and  time  neces¬ 
sary  to  achieve  a  desired  accuracy  level.  Traditional  finite  volume  methods  rely  on  extended 
stencils  to  achieve  high-order  accuracy,  which  leads  to  difficulty  in  achieving  stable  iterative 
algorithms  and  higher-order  accuracy  on  unstructured  meshes.  Alternatively,  an  attractive 
approach  for  achieving  higher-order  accuracy  is  the  discontinuous  Galerkin  (DG)  formu¬ 
lation  in  which  element-to-element  coupling  exists  only  through  the  fluxes  at  the  shared 
boundaries  between  elements. 

Recently,  Fidkowski  and  Darmofal  [22,  23]  developed  a  p-multigrid  method  for  the  solu¬ 
tion  of  high-order,  DG  discretizations  of  the  Euler  equations  of  gas  dynamics.  They  achieved 
significant  reductions  in  the  computational  time  required  to  obtain  high  accuracy  solutions 
by  using  high-order  discretizations  rather  than  highly  refined  meshes.  This  thesis  describes 
the  extension  of  the  algorithm  introduced  by  Fidkowski  and  Darmofal  to  viscous  flows. 


1.2  Background 

1.2.1  Higher-Order  Methods 

The  first  high-order  accurate  numerical  methods  were  spectral  methods  [24,  15],  where  the 
solution  of  a  differential  equation  is  approximated  over  the  entire  domain  using  a  high-order 
expansion.  Choosing  the  expansion  functions  properly,  one  can  achieve  arbitrarily  high- 
order  accuracy.  However,  because  of  the  global  nature  of  the  expansion  functions,  spectral 
methods  are  typically  limited  to  very  simple  domains  with  simple  boundary  conditions. 

Motivated  by  the  prospect  of  obtaining  the  rapid  convergence  rates  of  spectral  methods 
with  the  greater  geometric  versatility  provided  by  finite  element  methods,  researchers  in 
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the  early  1980s  introduced  the  p-type  finite  element  method.  In  the  p-type  finite  element 
method,  the  grid  spacing,  h,  is  fixed,  and  the  interpolation  order,  p,  is  increased  to  drive 
the  error  down.  In  1981,  Babuska  et  al.  [6]  applied  this  method  to  elasticity  problems. 
They  concluded  that  based  on  degrees  of  freedom,  the  rate  of  convergence  of  the  p-type 
method  cannot  be  slower  than  that  of  the  /i-type  and  that,  in  cases  with  singularities 
present  at  vertices,  the  convergence  rate  of  the  p- type  is  twice  as  fast.  In  1984,  Patera  [40] 
introduced  a  variant  p- type  method,  known  as  the  spectral  element  method,  and  used  it  to 
solve  the  incompressible  Navier-Stokes  equations  for  flow  in  an  expanding  duct.  Korczak 
and  Patera  [30]  later  extended  this  method  to  more  general,  curved  geometries. 

Finite  element  methods  are  attractive  for  achieving  high-order  accuracy  because,  for 
smooth  problems,  the  order  of  accuracy  is  controlled  by  the  order  of  the  solution  and  test 
function  spaces.  However,  it  is  well  known  that  standard,  continuous  Galerkin  methods  are 
unstable  for  the  convection  operator  [55].  Thus,  the  solution  of  the  Euler  or  Navier-Stokes 
equations  requires  the  addition  of  a  stabilization  term,  like  that  used  in  the  Streamwise 
Upwind  Petrov  Galerkin  method  [27]. 

Also  motivated  by  the  possibility  of  obtaining  spectral-like  results  in  a  more  flexible  geo¬ 
metric  framework,  Lele  [33]  introduced  up  to  tenth-order,  compact,  finite  difference  schemes. 
Using  a  Fourier  analysis  of  the  differencing  errors,  he  showed  that  the  compact,  high-order 
schemes  have  a  larger  resolving  efficiency,  where  resolving  efficiency  is  the  fraction  of  waves 
that  are  resolved  to  a  given  accuracy,  than  traditional  finite  difference  approximations.  This 
work  was  extended  to  more  general  geometries  by  Visbal  and  Gaitonde  [56],  who  used  up  to 
sixth-order,  compact  finite  difference  discretizations  to  solve  the  compressible  Navier-Stokes 
equations  on  curvilinear  meshes. 

Further  work  for  the  compressible  Navier-Stokes  equations  on  structured  meshes  is  pre¬ 
sented  by  Zingg  et  al.  [59],  who  compared  the  results  of  a  fourth-order  central  difference 
discretization  to  a  number  of  lower-order  schemes.  One  of  their  results  is  shown  in  Figure  1- 
1.  The  higher-order  scheme  is  seen  to  dramatically  outperform  the  second-order  (Matrix, 
CUSP,  and  Roe)  methods  in  terms  of  number  of  nodes  required  to  accurately  compute  the 
drag  on  an  airfoil  in  subsonic  flow. 

In  an  unstructured,  finite  volume  context,  Barth  [7]  introduced  the  fc-ex act  reconstruc¬ 
tion  method,  which  is  based  on  a  least-squares  reconstruction  procedure  that  requires  an 
extended  stencil.  As  noted  above,  the  large  stencil  typically  required  is  a  limiting  factor 
in  the  development  and  application  of  higher-order  finite  volume  methods.  Alternatively, 
Wang  et  al.  [57]  has  recently  developed  the  spectral  volume  method,  where  each  cell  (spectral 
volume)  in  the  domain  is  subdivided  into  additional  control  volumes.  The  state  averages 
within  these  control  volumes  are  used  to  build  a  higher-order  reconstruction  within  the 
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Figure  1-1:  Drag  results  for  subsonic  NACA  0012  test  case.  Taken  from  Zingg  et  al.  [59]. 
Reproduced  with  permission. 

spectral  volume.  To  increase  the  order  of  accuracy,  additional  degrees  of  freedom  are  added 
by  further  subdividing  each  spectral  volume.  Thus,  at  the  spectral  volume  level,  the  scheme 
has  a  nearest  neighbor  stencil  regardless  of  the  order  of  accuracy. 

1.2.2  Discontinuous  Galerkin  Methods 

In  1973,  Reed  and  Hill  [44]  introduced  the  DG  method  for  the  neutron  transport  equation. 
Since  that  time,  development  of  the  method  has  proceeded  rapidly.  Cockburn  et  al.  present 
an  extensive  history  of  DG  methods  in  [16].  Highlights  of  this  history  are  mentioned  here. 

In  1974,  LeSaint  and  Raviart  [34]  derived  the  first  a  priori  error  estimates  of  the  DG 
method  for  linear  hyperbolic  problems.  They  proved  a  rate  of  convergence  of  0{hP )  in 
the  L2(^)-norm-  Johnson  and  Pitkaranta  [29]  and  Richter  [45]  later  improved  upon  these 
original  estimates.  Johnson  and  Pitkaranta  proved  that,  in  the  most  general  case,  0(hP+ 1/2) 
is  the  optimal  convergence  rate,  while  Richter  showed  that,  assuming  the  characteristic 
direction  is  not  exactly  aligned  with  the  grid,  0(hP+1 )  can  be  obtained. 

A  breakthrough  in  the  application  of  DG  methods  to  nonlinear  hyperbolic  problems 
was  made  by  Cockburn  and  Shu  [18],  who  introduced  the  Runge  Kutta  Discontinuous 
Galerkin  (RKDG)  method.  The  original  RKDG  method  uses  an  explicit  TVD  second- 
order  Runge  Kutta  scheme  introduced  by  Shu  and  Osher  [49].  In  1989,  Cockburn  and 
Shu  [17]  generalized  the  method  to  be  higher-order  in  time  as  well  as  space.  Cockburn  and 
Shu  summarize  the  RKDG  method,  including  the  details  of  a  generalized  slope  limiter  for 
controlling  oscillations,  in  [20]. 

Independent  of  the  above  work,  Allmaras  [1]  and  Allmaras  and  Giles  [3]  developed  a 
second-order  DG  scheme  for  the  2-D  Euler  equations.  Their  method  is  the  extension  of 
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van  Leer’s  method  of  moments  [50]  from  the  1-D,  linear  wave  equation  to  the  2-D  Euler 
equations.  Thus,  it  requires  that  state  and  gradient  averages  be  computed  at  each  cell  to 
allow  linear  reconstruction  of  the  state  variables.  Halt  [25]  later  extended  this  work  to  be 
higher-order  accurate. 

For  elliptic  operators,  in  the  late  1970s  and  early  1980s,  Arnold  [5]  and  Wheeler  [58] 
introduced  discontinuous  finite  element  methods  known  as  penalty  methods.  While  these 
schemes  were  not  developed  as  DG  methods,  they  have  now  been  brought  into  the  unified 
DG  framework  [21].  More  recently,  many  researchers  [13,  9,  8,  20,  19,  42]  have  applied 
DG  methods  to  diffusive  operators.  One  procedure,  pioneered  by  Bassi  and  Rebay  [9,  11] 
and  generalized  by  Cockburn  and  Shu  [20,  19],  is  to  rewrite  a  second-order  equation  as 
a  first-order  system  and  then  discretize  the  first-order  system  using  the  DG  formulation. 
This  method  has  been  successfully  applied  to  the  compressible  Navier-Stokes  equations  and 
Reynolds  Averaged  Navier-Stokes  equations  by  Bassi  and  Rebay  [11,  10].  Arnold  et  al.  [21] 
provides  a  unified  analysis,  including  error  estimates,  of  most  of  the  DG  methods  available 
for  elliptic  operators. 


1.2.3  Multigrid  for  Aerodynamic  Applications 

The  use  of  multigrid  for  the  solution  of  the  Euler  equations  was  pioneered  by  Jameson  in 
1983,  who  demonstrated  a  significant  convergence  speedup  in  the  solution  of  2-D  transonic 
flows  on  structured  meshes  [28].  Since  that  time,  there  have  been  many  advances  in  the  ap¬ 
plication  of  multigrid  methods  to  aerodynamic  problems.  For  example,  Mavriplis  [36]  intro¬ 
duced  a  method  of  performing  multigrid  on  unstructured  triangular  meshes;  Allmaras  [2]  has 
examined  the  requirements  for  the  elimination  of  all  error  modes;  and  Pierce  and  Giles  [41] 
presented  efficient  methods  for  the  Euler  and  Navier-Stokes  equations.  Specifically  for  DG 
discretizations,  in  2002,  Bassi  and  Rebay  [12]  introduced  a  semi-implicit  p-multigrid  algo¬ 
rithm  and  used  it  to  solve  the  DG  discretization  of  the  Euler  equations. 

This  work  builds  directly  on  that  of  Fidkowski  and  Darmofal  [22,  23],  who  developed  a 
p-multigrid  algorithm  with  element-line  Jacobi  relaxation  for  solving  the  DG  discretization 
of  the  Euler  equations.  They  were  able  to  achieve  p-independent  asymptotic  convergence 
rates  as  well  as  significant  time  savings  versus  a  simpler  element  Jacobi  preconditioned  p- 
multigrid  scheme.  The  multigrid  solver  implemented  by  Fidkowski  is  used  here  with  only 
minor  modification  to  the  element  coupling  criterion  used  by  the  line  creation  procedure. 
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1.3  Outline  of  Thesis 


This  thesis  presents  a  multigrid  solution  technique  for  a  high-order  DG  discretization  of 
the  compressible  Navier-Stokes  equations.  The  equations  are  discretized  using  the  second 
method  of  Bassi  and  Rebay  (BR2)  [11,  10],  which  is  described  in  detail  in  Chapter  2. 
Chapter  3  describes  the  multigrid  algorithm  developed  by  Fidkowski  and  Darmofal  [22,  23] 
and  its  extention  to  high  Reynolds  number  viscous  flows.  Stability  analysis,  presented  in 
Chapter  4,  shows  that  the  single-step  element  and  element-line  Jacobi  relaxation  schemes  are 
stable  independent  of  p  and  flow  conditions.  Finally,  2-D  laminar  results  shown  in  Chapter  5 
demonstrate  that  higher-order  schemes  provide  significant  savings  in  terms  of  the  number 
of  elements,  degrees  of  freedom,  and  time  required  to  achieve  a  desired  accuracy  level. 
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Chapter  2 


Discontinuous  Galerkin 
Discretization 


This  chapter  develops  a  high-order  accurate  discretization  of  the  compressible  Navier-Stokes 
equations.  The  discretization  of  the  inviscid  terms  uses  the  standard  DG  formulation,  which 
relies  on  Riemann  solvers  for  the  calculation  of  the  inter-element  fluxes.  The  bulk  of  the 
chapter  focuses  on  the  discretization  of  the  viscous  terms,  which  is  done  using  the  second 
formulation  of  Bassi  and  Rebay  (BR2)  [11,  10]. 

2.1  DG  for  Euler 

The  two-dimensional  Euler  equations  in  strong,  conservation  form  are  given  by 


ut  +  V  •  Ti( u)  =  0, 


(2.1) 


where  u  is  the  conservative  state  vector, 

u  =  (p  pu  pv  pE)T , 
Ti  =  (F®,  Ff)  is  the  inviscid  flux  vector, 


(  pu  ^ 

^  pv  ^ 

Ff  = 

pu 2  +  p 

,  Ff  = 

puv 

% 

puv 

’  l 

pv2  +  p 

\  puH  ) 

\  pvH  y 
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p  is  the  fluid  density,  u  and  v  are  velocity  components,  p  is  the  pressure,  and  E  is  the 
total  internal  energy  per  unit  mass.  Thus,  the  total  enthalpy  per  unit  mass,  H ,  is  given  by 
H  =  E  +  p/ p,  and,  assuming  the  fluid  obeys  the  perfect  gas  equation  of  state,  the  pressure 
is  p  =  (7  —  1  )p[E  —  ( u 2  +  u2)/2],  where  7  is  the  ratio  of  specific  heats  of  the  fluid. 

Multiplying  Eqn.  2.1  by  a  vector- valued  test  function  v  and  integrating  by  parts,  one 
obtains  the  weak  formulation: 

f  v7u fdx—  f  VvT  f  vTEi  ■  n  ds  =  0,  Vv  G  i/1( 12), 

Jo.  Jfl  J  80. 

where  12  is  the  domain,  d  12  is  its  boundary,  and  n  is  the  outward  pointing  unit  normal.  To 
discretize  in  space,  define  Vp  to  be  the  space  of  discontinuous  vector-valued  polynomials 
of  degree  p  on  a  subdivision  Th  of  the  domain  into  non-overlapping  elements  such  that 
12  =  Uh-eTh  K-  Thus,  the  solution  and  test  function  space  is  defined  by 

K  =  {ve  L2(! 2)  I  v|*  €  Pp ,  Vk  €  Th},  (2.2) 


where  Pp  is  the  space  of  polynomial  functions  of  degree  at  most  p.  The  discrete  problem 
then  takes  the  following  form:  find  u/v  G  Vp  such  that  Vv/j  G  Vp, 


{  [  vh(uh)tdx-  f  Vv^  •  Pt  clx 

—  rj-\  ’J  K,  J  K, 

h 

+  [  v+THi(u+,u-,n)ds+  I 

J  dK,\dfl  J  d 


V+Tnb 

vh  '  Li 


(u+,utn)ds}  =  0,  (2.3) 


where  and  7i\{ u^,u^,n)  are  numerical  flux  functions  defined  for  interior 

and  boundary  faces,  respectively.  In  this  work,  the  Roe-averaged  flux  [46]  is  used  for  the 
inviscid  numerical  flux  on  interior  faces.  The  boundary  conditions  are  imposed  weakly  by 
constructing  an  exterior  boundary  state,  u£,  which  is  a  function  of  the  interior  state  and 
known  boundary  data.  Boundary  conditions  are  discussed  in  more  detail  in  Section  2.4.2. 
Furthermore,  the  (•)“*“  and  (•)“  notation  is  used  to  indicate  the  trace  value  taken  from  the 
interior  and  exterior  of  the  element,  respectively. 


2.2  DG  for  Navier-Stokes 

The  compressible,  two-dimensional  Navier-Stokes  equations  in  strong,  conservation  form 
are 

ut  +  V  •  Pi(u)  -  V  •  Pv(u,  Vu)  =  0,  (2.4) 
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where  the  conservative  state,  u,  and  inviscid  flux  vector,  J7,,  are  defined  in  Section  2.1.  The 
viscous  flux,  Tv  =  AVu  =  (FJ,  Fy),  is  given  by 


F 


X 

V 


Yy 


( 


V 

( 


dv 

dy 


§M2 1  -  %) 

Mlg  +  fJ) 

|)«+M|  +  ii)»+«s  / 

0  ' 

(du  _i —  9l>\ 

2  /nft)  _  9«\ 

S^dy  dx> 

g)’'  +  M(ft  +  g)»  +  « f  ) 


where  n  is  the  dynamic  viscosity  and  /t  is  the  thermal  conductivity. 


2.2.1  Flux  Formulation 

The  first  step  in  the  flux  formulation  is  to  rewrite  Eqn.  2.4  as  a  first-order  system.  To 
accomplish  this  reduction  of  order,  a  new  variable,  Q  =  A(,Vu,  is  defined.  Thus,  Eqn.  2.4 
can  be  written  as  a  first-order  system  in  terms  of  Q: 


U(  +  V  •  Ti  -  V  ■  Q  =  0, 


Q  -  AiVu  =  0. 


(2.5) 


Multiplying  Eqns.  2.5  by  test  functions  v  and  r,  respectively,  and  integrating  by  parts  gives 
the  weak  formulation: 

/  vru +  /  vT A  ■  nds  —  /  Vv7  •  A  dx 

Jo.  J  do.  J  n 

-  f  vTQ-hds+  f  VvT-Qdx  =  0,  VveE1^), 

Jan  Jn 

[  tt  ■  Qdx  -  f  uT(A^r)-nds 
Jn  JdQ 

+  f  urV  •  (Ay' t)  dx  =  0,  Vr  £  [iE1(fl)]2. 

A 

Then,  using  the  space  and  triangularization  Th  defined  in  Section  2.1,  the  spatial  dis¬ 
cretization  of  Eqn.  2.4  takes  the  following  form:  find  u^  £  and  Qh  £  [V^]2  such  that 
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Vvfc  £  Vph  and  Vrft  £  [V)(]2, 


kgth 


+ 


+ 


'  dnC\d£l 


'  dandfl 


v+THv{Ql,Qh)-nds 


El/ vh(uh)tdx-  [  Vv^-JTi0?x+  /  v+i'Hi(u+,ufe,n)ds 

_  nn  ^  J  Hi  J  hi  J  dn\d£l 

-  E  {  / 

J  KGTft  L  •'s'Al9n 

2h)  •  n^-s  -  /  Vv£  •  Qhdx}  =  0, 

E{  [ rl-Qhdx-  I  hu(u+,u~)T(A^Th)+  -hds 

^  J  n  J dn\d£l 

uhT(^-vrh)+  •  nds  +  /  u£V  •  t/0  rfx}  =  0, 

J  Hi  ' 


K&h 


'  duddfl 


(2.6) 


where  hu  and  are  numerical  fluxes  approximating  u,  and  HLV  and  H\,  are  numerical  fluxes 
approximating  ^Vu  on  interior  and  boundary  faces,  respectively.  Thus,  given  definitions 
of  these  numerical  fluxes,  Eqn.  2.6  gives  the  semi-discrete  flux  form  of  the  Navier-Stokes 
equations. 

Note  that  the  first  summation  over  k  in  Eqn.  2.6  is  exactly  the  same  as  that  in  the 
Euler  discretization  given  by  Eqn.  2.3.  These  terms  are  not  modified  by  the  addition  of  the 
viscous  discretization;  thus,  from  this  point  forward,  they  are  denoted  simply  as  E. 


2.2.2  Primal  Formulation 

To  simplify  the  notation  of  the  primal  form,  jump,  [•]],  and  average,  {•},  operators  are 
defined  for  interior  faces.  For  spatially  scalar  quantities,  the  operators  are  given  by 

[s]  =  s+n+ +  s_n“, 

W  =  M  +  *-), 

where  (-)+  and  (•)“  refer  to  trace  values  taken  from  opposite  sides  of  the  face.  Note  that 
the  unit  normals  n+  and  n~  are  outward  pointing  with  reference  to  the  sides  (-)+  and  (•)“, 
and,  thus,  n+  =  — n“.  Furthermore,  for  spatially  vector  quantities, 

M  =  <M-n+  +  M-n“, 

M  =  \iv+  +  v~)- 
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Via  substitution,  one  can  show  that, 


£  [  s+Tip+  ■  hds  =  [  [s|T  •  {</?}  ds  +  [  {s}T[v?]]ds,  (2.7) 

J  a  n\an  Jr.i  JTi 

where  F,;  is  the  union  of  all  interior  faces.  Applying  Eqn.  2.7,  Eqns.  2.6  become 


Vv£  •  Qh  dx 


“  /  brhY-{'Hv}ds 
dr; 


E+  E  [ 

Th  JK 

~  [  {v/i}T[Al  ds  -  f  VhTHbv  ■  hds  =  0, 

JTi  Jan 

22  [  j  Th  ‘  Qh  dyL  +  j  Uh  V  '  vTh )  dx  -  [h.a]T  •  {A^Th}  ds 

-  [  {hu}TlA^Thj  ds  -  [  u bhT(A^Th)+  •  n  ds  =  0. 

Jri  Jan 


k6  Th 


(2.8) 


Integrating  by  parts  and  using  Eqn.  2.7,  the  term  ^2KeTh  fK  uhT^  ■  (A^Th)  dx  in  Eqn  2.8 
can  be  rewritten: 

^22  J  uh^  ■  (AvTh)dx  =  (\uh\T -{A^t h}  +  {nh}T\A^ Th^jds 


kG  Th 


+  u+T(AvTh)+  ■  hds 

Jan 

22  I  rh'  (AVu/0  dx. 

_  s-rT1  J  /“v 


(2.9) 


KG  Th 


Then,  substituting  Eqn.  2.9  into  the  second  of  Eqns.  2.8  yields 


kG  Th 


22  {  I  rh-Qhdx-  j  rl  ■  (AVuk)  dx} 

^  }[hu  -  u,J]T  •  {A^rh}  +  {hu  -  u h}TlA^Th2j  ds 


~  /  (uh.  -  U-t)T{AvTh)+  ■  hds  =  0. 

Jan 


Finally,  defining  lifting  operators  S  and  Sl  as 


22  I  Th-Sdx  =  -  f  [hu-u,  h\T  ■  {A^Th}  ds  -  f  (ubh-u+)T{A^Th)+  -hds, 

J k  JTi  Jan 


k &ThJ  K  JTi 

22  [Th-dldx  =  -  f  {hu-uh}TlA^Thjds, 

k£Th  dTi 
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one  can  write  Qh  in  terms  of  u/n  <5,  and  8l: 


Qh  =  A.Vuft  -  5  -  <5*. 

Thus,  taking  t/j  =  Vv/,,  the  primal  form  is  given  by  the  following:  find  u E  Vj’  such  that 
Vvfe  E  Vph, 


E  +  ^2  J  Vv[  •  (AvVuh)dx  + J  -  uhjT  •  {^Vv/j}  -  [vft]]T  •  {Hv}}  ds 

K(E.Th  K  * 

^{hu  -  ufe}T[^VvftJ  -  {v4r[?it,]])  ds 


+ 


+ 


Ti 


Ian 


( uh  ~  uft)T(^VV/l)+  •  n  -  v+THtb  •  n 


ds  =  0. 


(2.10) 


Again,  to  fully  define  the  discretization,  one  must  define  the  numerical  fluxes  h?i(u4  u^j 
and  TLv(QJfl,  Q for  interior  faces  and  and  Ttbv  for  boundary  faces.  These  definitions  are 
considered  in  Section  2.2.3. 


2.2.3  Bassi  and  Rebay  Discretization 

Motivated  by  the  lack  of  any  upwinding  mechanism  in  the  viscous  terms  of  the  Navier-Stokes 
equations,  one  might  consider  using  central  fluxes  for  h„  and  7iv: 

h«  =  {u/J, 

nv  =  {Qh}. 

These  fluxes  were  originally  developed  and  applied  by  Bassi  and  Rebay  [8]  with  moderate 
success.  However,  they  are  used  here  only  to  motivate  the  final  flux  choice. 

Clearly,  both  fluxes  are  single-valued  on  a  given  face,  thus, 

{nv}  =  nv, 

iHvl  =  0, 

and 

=  o, 

=  -[Uhl- 


{h.u  -  uh } 
[hu  -  uh} 
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To  complete  the  interior  face  flux  definition,  note  that 


Hv  =  {Qh}  =  {AvVuh  -  6  -  S1}. 


However, 


uh}TlA^Thlds  =  0. 


Thus, 


Hv  =  {AVu/j}  -  {<5}. 


Finally,  boundary  conditions  are  of  the  form 


uh  =  u^(u+,BC  Data), 

U\  =  Hbv{ u+  Vu+,  BC  Data)  =  (A.VU/l)b  -  Sb. 

Thus,  the  final  form  of  the  discretization  is  as  follows:  find  u/,.  £  such  that  Vv^  £  Vf , 

E  +  ^2  [  Vvfe  '  (AVua)  dx  -  f  (K]t  ■  {^VvA}  +  [vh]T  •  {AVujj  ds 

k&Th  '  Ti 

+  [  [K]T  •  {<5}  ds  +  [  (uh  -  u£)T(^Vvh)+  ■  hds 

JVi  Jan 

—  [  v^T(AvV\ih)b  ■  hds  +  f  v^T Sb  ■  hds  =  0. 

Jan  Jan 

Unfortunately,  as  shown  by  numerous  authors  [20,  21,  42],  this  discretization  is  problem¬ 
atic  for  multiple  reasons.  First,  on  some  meshes  the  viscous  contribution  to  the  Jacobian 
may  be  singular.  While  the  inviscid  terms  should  make  the  discrete  system  non-singular, 
this  result  is  clearly  undesirable.  Second,  stability  and  optimal  order  of  accuracy  in  L2 
cannot  be  proven  [21].  In  fact,  Cockburn  and  Shu  [20]  have  shown  that  for  purely  elliptic 
problems,  odd  order  interpolants  produce  sub-optimal  order  of  accuracy  equal  to  0(hp). 
Finally,  the  scheme  is  not  compact.  This  non-compactness  is  introduced  through  the  global 
lifting  operator  <5,  as  illustrated  in  1-D  by  Figure  2-1.  The  figure  shows  that,  while  the  resid¬ 
ual  on  element  k,  Rk,  depends  explicitly  on  only  and  5  on  the  element  and  its  neighbors 
(elements  k  —  1,  k,  and  k  +  1),  the  global  lifting  operators  on  the  neighboring  elements  also 
depend  on  their  neighbors.  Thus,  the  dependence  of  Rk  is  extended  to  elements  k  —  2  and 
k  +  2  if  the  additional  degrees  of  freedom  for  8  are  eliminated. 

Given  the  drawbacks  of  the  previous  scheme,  a  modification  of  this  scheme  which  will 
make  it  stable,  optimally  convergent,  and  compact  is  desired.  Bassi  and  Rebay  [11,  10]  have 
proposed  modified  numerical  fluxes  which  accomplish  these  goals  by  replacing  the  global 
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Figure  2-1:  1-D  stencil  for  first  Bassi  and  Rebay  scheme 


lifting  operator  with  a  local  lifting  operator.  The  local  lifting  operator,  or  auxiliary  variable, 
6f,  is  defined  by  the  following  problem:  find  <5/  E  [Vj[]2  such  that  Vr^  E  [V((]2, 

[  T%-6fdx=f  [u hjT-{A^Th}ds  (2.11) 

J  J  <J  f 

for  interior  faces,  and 

[  Th  ■  <5/dx  =  -  [  (ubh-  u+)T[(A^Th)  •  n]+  da,  (2.12) 

J  K+  J  CTj 

for  boundary  faces,  where  cry  denotes  a  single  interior  face  and  crb  denotes  a  single  boundary 
face. 

Replacing  the  global  lifting  operator,  S,  with  the  local  lifting  operator,  Sf,  and  mul¬ 
tiplying  by  a  stabilization  parameter,  rjf,  which  is  discussed  in  Section  2.3,  the  numerical 
fluxes  TLV  and  Hb  become 


TLV  =  {A.Vuft}  -  r)f{6f}, 

Hbv  =  (AvVuh)b -VfSb. 

The  hM  and  u bh  numerical  fluxes  are  not  modihed.  Thus,  substituting  the  new  numerical 
fluxes  into  Eqn.  2.10  gives  the  second  Bassi  and  Rebay  discretization:  find  u/j  E  such 
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Figure  2-2:  1-D  stencil  for  BR2  scheme 


that  Vv/j  G  Vf  , 


E  +  ^2  [  Vvfe  '  (AVuft)  dx  -  f  (juJT  •  {iJ’Vvh}  +  [vh]T  •  {i„Vu,,}j  ds 


Ti 

,b  „  +  \T /  A Tt 


+  /  [vj  •  77/{<5/}ds+  /  (u  -u^n^Vv^+mds 

Jr;  Van 

-  [  v^T(AvVuh)b  ■  hds  +  [  v^Tiifdbf  ■  iids  =  0. 

Jan  Jan 


(2.13) 


Note  that,  unlike  the  original  choice  of  fluxes,  the  BR2  form  has  a  compact  stencil,  meaning 
that  only  elements  that  share  a  face  are  coupled.  Figure  2-2  shows  the  stencil  of  the  BR2 
scheme  in  1-D.  It  shows  that  the  compact  stencil  follows  from  the  fact  that  the  local  lifting 
operators,  Sf,  at  a  given  face  depend  only  on  the  elements  that  share  that  face. 

Furthermore,  for  purely  elliptic  problems  with  homogeneous  Dirichlet  boundary  condi¬ 
tions,  Arnold  et  al.  [21]  proves  optimal  error  convergence  in  the  L 2  norm  for  p  >  1  when 
rjf  >  3  (note:  the  condition  rjf  >  3  is  required  to  prove  stability).  However,  in  Section  2.3 
we  propose  a  definition  of  rjf  which  is  generally  less  than  this  value  but  is  required  to  pro¬ 
duce  optimal  accuracy  for  p  =  0;  in  practice,  we  have  not  found  that  this  lower  value  of  i]f 
causes  a  loss  of  stability. 

While  the  BR2  discretization  is  certainly  not  the  only  scheme  proposed  in  the  literature 
which  achieves  optimal  order  of  accuracy,  it  has  some  clear  advantages.  Other  schemes 
that  achieve  optimal  order  of  accuracy  include  the  local  discontinuous  Galerkin  (LDG) 
method  [19,  20],  a  penalty  method  proposed  by  Brezzi  [42],  and  the  Baumann  and  Oden 
scheme  [13].  However,  each  of  these  has  some  drawbacks.  LDG  is  not  compact  on  general, 
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unstructured  meshes.  The  stabilization  parameter  required  for  optimal  accuracy  in  the 
Brezzi  scheme  can  grow  very  large  at  high  p  because  pf  ~  h~2p ,  and  the  Baumann  and 
Oden  scheme  is  only  stable  for  p  >  2.  In  fact,  the  BR2  scheme  is  the  only  one  proposed  in 
the  literature  to  achieve  optimal  order  of  accuracy  for  all  p  >  1  with  a  compact  stencil. 

The  final  discrete  form  is  completed  by  choosing  a  basis  and  time  discretization,  as 
shown  in  Section  2.5. 

2.3  The  Stabilization  Parameter,  ijf 

As  discussed  in  Section  2.2.3,  Pf  >  3  is  required  to  prove  stability  for  the  BR2  scheme.  For 
p  >  1,  there  are  no  other  conditions  on  the  choice  of  pf.  However,  for  p  =  0,  if  pj  is  not 
set  appropriately,  the  error  may  not  converge  with  h.  A  convergent  p  =  0  discretization  is 
desired  such  that  it  may  be  used  as  the  coarsest  order  in  the  p-multigrid  algorithm.  This 
section  shows  the  derivation  of  the  choice  of  pj  used  in  this  work. 

2.3.1  Motivation 

To  motivate  the  derivation  of  a  non-constant  p f ,  results  from  the  solution  of  Poisson’s 
equation  using  Pf  =  3  are  shown.  Thus,  the  problem  is 

-V-(Vu)  =  /  in  Q  (2.14) 

u  =  g  on  dll, 

where 

/  =  -6(x  +  y), 

g  =  x3  +  y3. 

For  brevity,  the  BR2  discretization  of  this  problem  is  not  given  here.  However,  it  is  given 
for  g  =  0  in  Eqn.  2.15. 

To  conduct  a  mesh  refinement  study,  the  domain,  a  unit  square  whose  lower  left  corner 
lies  on  the  origin,  is  triangulated  into  grids  of  44,  176,  and  704  elements.  Results  of  the  study 
are  shown  in  Figure  2-3.  For  p  >  1,  the  L2-nornr  of  the  solution  error  converges  optimally 
at  a  rate  of  0{HP+1).  However,  for  p  =  0,  the  L2-nornr  of  the  error  is  approximately 
constant  with  grid  refinement.  This  behavior  results  from  the  fact  that  for  p  =  0,  only 
the  auxiliary  variable  terms  remain  on  the  left  hand  side  of  the  BR2  form.  If  pf  is  not  set 
appropriately,  the  p  =  0  discretization  is  not  consistent.  Thus,  a  definition  of  pf  that  makes 
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Figure  2-3:  Mesh  refinement  results  for  0  <  p  <  2  using  rjf  =  3 


the  discretization  consistent  is  sought. 


2.3.2  Formulation 


For  simplicity,  the  analysis  is  shown  for  Poisson’s  equation  with  homogeneous  Dirichlet 
boundary  conditions.  Thus,  the  problem  is  given  by  Eqn.  2.14  with  g  =  0.  For  this  case, 
the  BR2  discretization  is  given  by  Eqn.  2.15:  find  Uh  E  Vj'  such  that  Vvh  E  V^, 


F,  Vu/j  •  Vr^dx 

K&Th 


(M  •  {Vu/j}  +  [uj  •  (Vrt/j)  ds 


+  /  Vflvh}  •  {$f}ds+  /  -u^Vv^-hds 
d  Ti  Jon 


ion 


an 

u+Vu+-nds  +  f  r,fv+Sbrhds=  ^  [  vhf  dx, 

dsn  kgT^  d/C 


(2.15) 


where  is  now  the  space  of  discontinuous  scalar-valued  polynomials  of  degree  p  on  a 
subdivision  T),  of  the  domain  into  non-overlapping  elements  such  that  Id  = 

For  p  =  0,  Vuft  =  Vvh  =  0.  Thus,  Eqn.  2.15  becomes 


dfhhj- {Sf}ds+  r)fv£dbf -nds  = 


on 


E 

/ten 


L 


Vhf  dx. 


(2.16) 
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Examining  the  test  function  Vhi  associated  with  element  t Cj, 


l  if  x  e  m 
0  if  x  ^  Ki, 


Eqn.  2.16  reduces  to 


(2.17) 


where  Ki  has  no  faces  on  the  domain  boundary.  Using  the  definition  of  the  auxiliary  variable 
given  in  Eqn.  2.11,  a  p  =  0  basis  for  Tf l  yields 


(Srh)+  =  -uh)sa,  (2.18) 

(Srh)~  =  uhW,  (2.19) 


on  a  given  edge,  a  e  chcj,  of  length  sa.  A+  and  Aa  are  the  areas  of  the  elements  adjoining 
edge  a.  Substituting  Eqns.  2.18  and  2.19  into  Eqn.  2.17  gives 


E 

a^di^i  L 


^fsl(u+-uh) 


A+  +  An 


(2.20) 


To  motivate  the  definition  of  rjf  for  p  =  0,  apply  the  Divergence  theorem  to  Poisson’s 
equation, 

/  f  cht  =  —  V  ■  (Vu)  dx  =  —  /  Vu  •  hds. 

J  Ki  J  J  &Ki 


For  each  edge,  define  an  approximation  to  the  directional  derivative  in  the  n  direction  to 
be 


„  „  uh 

Vu  •  n  =  — 

where  An  is  the  distance  between  the  +  and 
direction.  For  triangles, 


An  ’ 

—  element  centroids  projected  onto  the  n 


=  \{K  +  K), 


where  h±  denotes  the  element  height  for  elements  adjacent  to  face  a,  as  shown  in  Figure  2-4. 
Thus, 


Vu  ■  hds  = 


crtzdKi 


.+ 


An 


(2.21) 
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Figure  2-4:  Definitions  of  h+,  h  ,  and  An 


Then,  substituting  Eqn.  2.21  into  Eqn.  2.20  and  solving  for  rjf  yields 

4  A+A- 
Vf  saAna  At  +  A~  ' 

Eqn.  2.22  is  valid  for  general  elements.  For  triangles,  A±  =  ^h±sa,  thus, 


(2.22) 


rjf  = 


i  _i_  i  ( Ai  i  Al  \ 
1  "T  2U-  +  h+) 


An  analogous  procedure  shows  that,  for  triangular  elements  with  straight  edges,  on  the 
boundary,  r)f  =  3/2  always. 

Thus,  using  this  definition,  r]f  <  3  for  every  face.  This  result  appears  problematic  given 
that  the  stability  proof  requires  r]f  >  3.  However,  in  practice,  no  problems  have  arisen. 
Results  for  the  case  examined  in  Section  2.3.1  are  given  in  Figure  2-5.  The  figure  shows 
that  the  L2  norm  of  the  solution  is  converging  at  0(hp+1)  for  0  <  p  <  2.  Furthermore,  for 
compressible  Navier- Stokes,  Bassi  and  Rebay  [11,  10]  successfully  applied  the  discretization 
with  rjf  =  1  for  p  >  1,  and  it  is  used  here  with  r]f  as  given  in  Eqn.  2.22  without  detrimental 
effects.  This  fact  is  shown  by  the  accuracy  studies  presented  in  Chapter  5. 


2.4  Boundary  Treatment 

2.4.1  Geometry  Representation:  Curved  Boundaries 

As  shown  by  Bassi  and  Rebay  [9],  high-order  DG  methods  are  highly  sensitive  to  the 
geometry  representation.  Thus,  it  is  necessary  to  build  a  higher-order  representation  of  the 
domain  boundary.  In  this  work,  the  geometry  is  represented  using  a  nodal  Lagrange  basis. 
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Figure  2-5:  Mesh  refinement  results  for  0  <  p  <  2  using  r]f  as  defined  in  Eqn  2.22 

Thus,  the  mapping  between  the  canonical  triangle  and  physical  space  is  given  by 

x  = 

j 

where  (j>j  is  the  jth  basis  function,  £  is  the  location  in  the  reference  space,  and  Xj  is  the 
location  of  the  jth  node  in  physical  space.  In  general,  the  Jacobian  of  this  mapping  need 
not  be  constant,  meaning  that  triangles  with  curved  edges  can  be  mapped  to  the  straight 
edged  canonical  element.  Thus,  by  placing  the  non-interior,  higher-order  nodes  on  the  real 
domain  boundary,  a  higher-order  geometry  representation  is  achieved. 

Two  notes  about  this  geometry  representation  must  be  made.  First,  on  a  curved  element, 
basis  functions  that  are  polynomials  of  order  p  on  the  canonical  element  are  not  polynomials 
of  order  p  in  physical  space.  Thus,  the  interpolation  order  in  physical  space  is  less  than  p. 
Second,  as  discussed  in  Section  5.2,  it  is  unclear  if  this  choice  of  geometry  representation 
is  optimal  as  oscillations  in  the  interpolated  geometry  may  have  detrimental  effects  on  the 
order  of  accuracy. 

2.4.2  Boundary  Conditions 

Boundary  conditions  are  enforced  weakly  via  the  domain  boundary  integrals  appearing  in 
Eqn.  2.13.  To  evaluate  these  integrals,  two  previously  undefined  terms  must  be  defined: 
and  (AvVuh)b-  This  section  defines  these  terms  for  various  boundary  conditions. 
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Full  State  Condition 


In  some  circumstances,  the  entire  state  vector  at  the  boundary,  u^,  may  be  known.  In  these 
cases,  the  inviscid  flux  is  computed  using  the  Riemann  solver  exactly  as  if  the  face  were  an 
interior  face: 

=  Wi(u£,ufc,n). 

No  conditions  are  set  on  the  viscous  flux,  thus,  (At,Vu/j)b  is  set  by  interpolating  (^0Vu(,)+ 
to  the  boundary,  Sb  is  computed  using  as  shown  in  Eqn.  2.12. 


Inflow/Outflow  Conditions 

At  an  inflow/outflow  boundary,  the  boundary  state,  u^,  is  defined  using  the  outgoing  Rie¬ 
mann  invariants  and  given  boundary  data.  Table  2.1  details  the  inflow/outflow  conditions 
used  in  this  work.  Note  that,  while  the  Euler  equations  are  well-posed  with  just  the  bound- 


Table  2.1:  Types  of  inflow/outflow  boundary  conditions. 


Condition 

Number  of  BCs 

Value  Specified 

Outgoing  Invariant 

Subsonic  Inflow 

3 

Tt,  PTi  ot 

J+ 

Supersonic  Inflow 

4 

P ,  pu,  pv,  pE 

None 

Subsonic  Outflow 

1 

V 

J+,  v,  s 

Supersonic  Outflow 

0 

None 

J+,  J~ ,  v,  s 

ary  conditions  listed  in  Table  2.1,  the  Navier-Stokes  equations  are  not.  Conditions — numeric 
and/or  physical — are  needed  to  set  (At,V Uh)b  and  Sb.  For  this  work,  (A„V u/,)fc  is  set  by 
extrapolating  from  the  interior  of  the  domain,  and  Sb  is  computed  as  shown  in  Eqn.  2.12. 
These  gradient  conditions  have  not  been  theoretically  investigated  and  may  be  expected  to 
degrade  accuracy  and  stability  at  low  Re. 

No  Slip  Wall  Conditions 

Two  no  slip  wall  conditions  are  used:  adiabatic  and  isothermal.  At  a  no  slip,  isothermal 
wall,  the  velocity  components  and  static  temperature  are  set.  Thus, 

ub  =  0, 
vb  =  0, 

Tb  rp 

—  -I- wall  • 
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To  set  the  full  state,  these  conditions  are  combined  with  the  static  pressure,  p,  which  is 

computed  using  the  interior  state  interpolated  to  the  boundary  and  given  boundary  data 
(„,b  _  n,b  _  nV 


\ 

No  physical  conditions  are  set  on  the  viscous  flux  at  the  wall.  However,  the  scheme  requires 
that  this  flux  be  set.  Thus,  it  is  extrapolated  from  the  interior.  The  auxiliary  variable,  Sf, 
is  computed  as  stated  in  Eqn.  2.12. 

At  a  no  slip,  adiabatic  wall,  the  velocity  components  and  the  heat  transfer  to  the  wall 
are  given: 


P=  (7-  1  )pE+. 

Thus,  the  boundary  state  is  given  by 


ub  =  0, 

vb  =  0, 
dT 

7T  =  °- 

dn  b 

Only  two  conditions  on  the  boundary  state  are  known,  thus,  two  variables  must  be  set  using 
interior  data  and  used  to  compute  the  boundary  state.  Static  pressure,  p.  and  stagnation 
enthalpy  per  unit  mass,  H ,  have  been  chosen: 


P 

H 


(7-1  )PE+, 
pE+  +  p 
P+ 


Thus,  the  boundary  state  is  given  by 


uh 


/  Pb  \ 

(  7  V  \ 

7-1  H 

pub 

0 

pvb 

0 

,  nPb  , 

,  nP+  , 

The  adiabatic  condition,  combined  with  the  no  slip  condition,  requires  the  viscous  flux 
associated  with  the  energy  equation — the  fourth  flux  component — to  be  zero.  The  other 
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components  of  the  viscous  flux  are  unconstrained  by  this  condition.  Therefore,  the  interior 
viscous  flux  is  used  for  these  components  and,  thus,  the  boundary  viscous  flux  is  given  by 


[(XV Uh)b  •  n] 


^  [(XVuft)+  •  n]r  ^ 
[(A,Vu/l)+  •  h]2 
[(A,Vu/l,)+  •  n]3 

V  0 


Furthermore,  since  the  energy  equation  viscous  flux  is  specified,  the  corresponding  compo¬ 
nent  of  the  local  lifting  operator  is  set  to  zero.  All  other  auxiliary  variables  are  computed 
in  the  usual  fashion.  Thus, 


sf 


•  n  = 


(  [<*/ -A]  i  N 
[Sbf  •  n]2 
[Sbf  •  n]3 

V  0 


2.5  Final  Discrete  System 

To  define  the  final  discrete  form,  it  is  necessary  to  select  a  basis  for  the  space  Vj'  and 
discretize  in  time.  For  the  basis,  a  set  of  element-wise  discontinuous  functions,  { cj)j },  is 
chosen  such  that  each  (f>j  has  local  support  on  one  element  only.  Thus,  the  discrete  solution 
has  the  form 

j 

In  this  work,  a  nodal  Lagrange  basis  with  uniformly  spaced  nodes  is  used.  However,  it 
should  be  noted  that  this  basis  becomes  poorly  conditioned  as  p  increases,  which  can  degrade 
the  iterative  convergence  rate.  A  hierarchical  basis  like  that  used  in  [22]  can  eliminate  this 
problem  with  the  additional  benefit  of  simplifying  the  multigrid  prolongation  and  restriction 
operators. 

A  backward  Euler  discretization  is  used  for  the  time  integration.  Thus,  the  discrete 
system  is  given  by 

-^M(un+1  -  un)  +  R(u”+1)  =  0, 

where  A4  is  the  mass  matrix  and  R(un+1)  is  the  steady  residual  vector. 

This  work  is  principally  concerned  with  steady  problems.  However,  the  unsteady  term 
is  included  to  improve  the  performance  of  the  solver  in  the  initial  iterations.  After  this 
initial  transient  period,  At  — >  oo  [22] . 
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Chapter  3 


Solution  Method 


Applying  the  discretization  developed  in  Chapter  2,  the  discretized  compressible  Navier- 
Stokes  equations  are  given  by  a  nonlinear  system  of  equations  R(u)  =  0.  To  solve  this 
system,  a  p-multigrid  scheme  with  element-line  Jacobi  smoothing,  developed  by  Fidkowski 
and  Darmofal  [22,  23],  is  used.  A  general  preconditioned  iterative  scheme  can  be  written 

un+i  =u^_p-lR(u^)i 

where  the  preconditioner  matrix,  P,  is  an  approximation  to  the  Jacobian,  Two  types  of 
preconditioners  are  examined:  element  Jacobi,  where  the  unknowns  on  a  single  element  are 
solved  simultaneously,  and  element-line  Jacobi,  where  the  unknowns  on  a  line  of  elements 
are  solved  simultaneously.  Details  of  both  smoothers  as  well  as  the  line  creation  procedure 
and  multigrid  solver  are  presented.  The  discussion  draws  heavily  on  [22]. 


3.1  Preconditioners 


For  the  element  Jacobi  scheme,  the  unknowns  on  a  single  element  are  solved  simultaneoulsy. 
The  diagonal  blocks  of  the  Jacobian  matrix  represent  the  influence  of  the  state  variables 
in  a  given  element  on  the  residual  in  that  element.  Thus,  the  preconditioner  matrix  is  the 
block  diagnonal  of  the  Jacobian.  To  improve  the  robustness  during  the  initial  iterations, 
the  block  diagonal  is  augmented  by  an  unsteady  term.  Thus,  the  preconditioner  is, 


p  =  d  +  a iM- 


where  D  is  the  block  diagnonal  of  ^  and  Ai  is  the  mass  matrix.  As  the  solution  converges, 
At  — ►  oo. 
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The  addition  of  the  unsteady  term  does  not  change  the  block  diagonal  structure  of  P. 
Thus,  it  is  inverted  one  block  at  a  time  using  Gaussian  elimination. 

The  element-line  Jacobi  scheme  is  slightly  more  complex.  In  strongly  convective  systems, 
transport  of  information  proceeds  along  characteristic  directions.  By  solving  implicitly  on 
lines  of  elements  connected  along  these  directions,  one  can  alleviate  the  stiffness  associated 
with  strong  convection.  Furthermore,  for  viscous  flows,  the  element-line  Jacobi  solver  is  an 
important  ingredient  in  removing  the  stiffness  associated  with  regions  of  high  grid  anisotropy 
frequently  required  in  viscous  layers  [2,  37].  Thus,  the  element-line  Jacobi  scheme  requires 
the  ability  to  construct  lines  of  elements  based  on  some  measure  of  element-to-element 
coupling  and  to  solve  implicitly  on  each  line.  This  section  considers  the  construction  and 
inversion  of  the  preconditioner  matrix  given  a  set  of  lines.  The  line  creation  algorithm  is 
described  in  Section  3.2. 

Given  a  set  of  IV;  lines,  the  preconditioner  matrix,  P,  is  composed  of  TV;  tridiagonal 
systems  constructed  from  the  linearized  flow  equations.  Let  the  tridiagonal  system  for  line 
l ,  where  1  <  l  <  IV;,  be  written  M;.  Then,  denoting  the  number  of  elements  in  line  l  as  n;, 
M;  is  a  block  rq  x  n;  matrix.  As  before,  the  on-diagonal  blocks  represent  the  influence  of  the 
state  variables  in  a  given  element  on  the  residual  in  that  element.  The  off-diagonal  blocks, 
M;(j,  k),  represent  the  influence  of  the  states  in  element  k  on  the  residual  in  element  j. 

As  in  the  element  smoother  case,  the  element-line  smoother  is  augmented  by  an  unsteady 
term  to  improve  robustness.  Thus,  the  final  form  of  the  preconditioner  is, 

where  M  is  the  set  of  assembled  M;  matrices. 

Inversion  of  P  uses  a  block-tridiagonal  algorithm  in  which  the  diagonal  block  is  LU 
decomposed.  As  the  dominant  cost  of  the  element-line  Jacobi  solver  (especially  for  higher- 
order  schemes)  is  the  LU  decomposition  of  the  diagonal,  the  computational  cost  of  the 
element-line  Jacobi  smoother  scales  as  that  of  the  simpler  element  Jacobi.  However,  the 
performance  of  the  element-line  Jacobi  smoother  is  significantly  better  due  to  the  increased 
implicitness  along  strongly  coupled  directions. 


3.2  Line  Creation 

The  effectiveness  of  the  element-line  Jacobi  smoother  depends  largely  on  the  quality  of  lines 
produced  by  the  line  creation  procedure.  For  inviscid  or  nearly  inviscid  flows,  information 
flows  along  characteristic  directions  set  by  convection.  Thus,  the  lines  should  be  aligned 
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with  the  streamline  direction  to  alleviate  the  stiffness  associated  with  strong  convection. 
In  viscous  flows,  the  effects  of  diffusion  and  regions  of  high  grid  anisotropy,  such  as  those 
typically  associated  with  boundary  layers,  couple  elements  in  directions  other  than  the 
convection  direction.  Connecting  elements  in  these  off-convection  directions  can  alleviate 
the  stiffness  associated  with  diffusion  and  grid  anisotropy. 

The  line  creation  procedure  is  divided  into  two  parts:  the  connectivity  criterion,  which  is 
a  measure  of  the  coupling  between  elements,  and  the  line  creation  algorithm,  which  connects 
elements  into  lines  based  on  the  connectivity  criterion. 


3.2.1  Connectivity  Criterion 

The  measure  of  coupling  used  in  this  work  is  similar  to  that  used  in  the  nodal  line  creation 
algorithm  of  Okusanya  [39].  In  that  algorithm,  the  coupling  was  taken  directly  from  the 
discretization. 

In  this  work,  the  coupling  is  based  on  a  p  =  0  discretization  of  the  scalar  transport 
equation, 

V  •  (pu(j>)  —  V  •  (/uV <f)  =  0, 


where  pu  and  p  are  taken  from  the  solution  at  the  current  iteration.  More  specifically,  the 
coupling  between  two  elements  j  and  k  that  share  a  face  is  given  by 


Jj,k 


=  max 


( 

dRj 

dRk 

V 

d(pk 

1 

d<t>j 

While  this  definition  of  coupling  does  not  represent  the  exact  coupling  between  ele¬ 
ments  for  the  higher-order  Navier-Stokes  discretization  being  solved,  it  captures  the  rele¬ 
vant  features — the  effects  of  convection  and  diffusion — and  remains  unique.  In  the  p  =  0 
discretization  of  the  scalar  transport  equation,  the  off-diagonal  components  of  the  Jacobian 
matrix  are  scalars.  For  a  p  >  1  discretization  of  the  scalar  transport  equation  or  any  order 
discretization  of  a  system  of  equations,  the  off-diagonal  blocks  of  the  Jacobian  are  matrices. 
Thus,  a  matrix  norm  would  be  required  to  make  the  coupling  unambiguous. 


3.2.2  Line  Creation  Algorithm 

After  computing  the  elemental  coupling,  lines  of  elements  are  formed  using  the  line  creation 
algorithm  developed  by  Fidkowski  and  Darnrofal  [22,  23].  The  summary  given  here  is  based 
on  [22], 

The  line  creation  process  is  divided  into  two  stages:  line  creation  and  line  connection. 
Let  N(j;  /)  denote  the  element  adjacent  to  element  j  across  face  /,  and  let  F(j)  denote  the 
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set  of  faces  enclosing  element  j.  Then,  the  line  creation  algorithm  is  as  follows: 

Stage  I:  Line  Creation 

1.  Obtain  a  seed  element  i 

2.  Call  MakePath(i)  -  Forward  Path 

3.  Call  MakePath(i)  -  Backward  Path 

4.  Return  to  (1).  The  algorithm  is  finished  when  no  more  seed  elements  exist. 

MakePath(j) 

While  path  not  terminated: 

For  element  j .  pick  the  face  /  E  F(j)  with  highest  connectivity,  such  that  element 
k  =  N(j ;  /)  is  not  part  of  the  current  line.  Terminate  the  path  if  any  of  the 
following  conditions  hold: 

-  face  /  is  a  boundary  face 

-  element  k  is  already  part  of  a  line 

-  C(j,  k )  is  not  one  of  the  top  two  connectivities  in  element  k 
Otherwise,  assign  element  j  to  the  current  line,  set  j  =  k,  and  continue. 

After  completing  Stage  I,  it  is  possible  that  the  endpoints  of  two  lines  are  adjacent  to  each 
other,  as  illustrated  in  Figure  3-1.  Elements  a  and  b  have  not  been  connected  because  the 


Figure  3-1:  Possible  line  configuration:  (a)  after  Stage  I  and  (b)  after  stage  II.  Reproduced 
with  permission  from  [22]. 
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connectivity,  Ca  j,,  is  the  minimum  connectivity  for  both  elements.  However,  the  pairs  a,  c 
and  b,  d  have  not  been  connected  because  Ca,c  and  are  the  minimum  connectivities  for 
elements  c  and  d,  respectively.  For  best  solver  performance,  it  is  desirable  to  use  lines  of 
maximum  length.  Thus,  it  is  necessary  to  connect  elements  a  and  b.  The  line  connection 
stage  accomplishes  this  goal. 

Stage  II:  Line  Connection 

1.  Loop  through  endpoint  elements,  j .  of  all  lines.  Denote  by  Hj  C  F(j)  the  set  of  faces 
h  of  j  that  are  boundary  faces  or  that  have  N(j\  h)  as  a  line  endpoint. 

2.  Choose  h  £  Hj  of  maximum  connectivity.  If  h  is  not  a  boundary  face,  let  k  =  N(j  \  h). 

3.  If  k  has  no  other  neighboring  endpoints  of  higher  connectivity,  and  no  boundary  faces 
of  higher  connectivity,  then  connect  the  two  lines  to  which  j  and  k  belong. 

Proofs  that  both  stages  of  the  line  creation  algorithm  result  in  a  unique  set  of  lines,  inde¬ 
pendent  of  seed  element,  are  provided  in  [22]. 

The  lines  formed  for  a  2-D  NACA  0012  test  case  are  shown  in  Figure  3-2.  As  shown 


(a)  Outer  flow  (b)  Trailing  edge 


Figure  3-2:  Lines  formed  around  NACA  0012  in  M  =  0.5,  Re  =  5000,  a  =  0°  flow 

in  Figure  3-2(a),  the  lines  in  the  outer  flow,  where  convection  dominates,  simply  follow  the 
streamline  direction.  In  viscous  regions — the  boundary  layer  and  wake  in  this  problem — the 
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effects  of  high  aspect  ratio  elements  and  diffusion  produce  elements  that  are  tightly  coupled 
in  the  direction  normal  to  convection.  Thus,  as  shown  in  Figure  3-2(b),  lines  are  formed 
normal  to  the  streamline  direction.  Details  and  results  for  this  case  are  given  in  Chapter  5. 


3.3  j9-Multigrid 

3.3.1  Motivation 

The  use  of  multigrid  techniques  is  motivated  by  the  observation  that  the  smoothers  devel¬ 
oped  in  Section  3.1  are  ineffective  at  eliminating  low-frequency  error  modes  on  the  fine  grid. 
In  standard  /i-multigrid,  spatially  coarser  grids  are  used  to  correct  the  error  on  the  fine  grid. 
On  a  coarser  grid,  the  low-frequency  error  modes  from  the  fine  grid  appear  as  high-frequency 
modes,  and,  thus,  the  smoothers  can  effectively  correct  them.  p-Multigrid  uses  the  same 
principle  except  that  lower-order  approximations  serve  as  the  “coarse  grid”  [48,  26]. 

Furthermore,  p-multigrid  fits  naturally  into  the  high-order  DG,  unstructured  grid  frame¬ 
work.  Unlike  /i-multigrid,  spatially  coarser  meshes  are  not  required.  Thus,  no  element 
agglomeration  or  re-meshing  procedures  are  necessary.  Only  prolongation  and  restriction 
between  orders  are  required.  Moreover,  the  prolongation  and  restriction  operators  are  local, 
meaning  that  they  must  only  be  stored  for  the  canonical  element. 

3.3.2  FAS  and  Two-level  Multigrid 

The  multigrid  method  used  here  is  the  Full  Approximation  Scheme  (FAS),  introduced  by 
Brandt  [14].  The  following  description  is  taken  from  Fidkowski  [22], 

Consider  the  discretized  system  of  equations  given  by 

RP(UP)  =  fP 

where  up  is  the  discrete  solution  vector  for  pth  order  interpolation  on  a  given  grid,  Rp(up) 
is  the  associated  nonlinear  system,  and  fp  is  a  source  term  (zero  for  the  fine-level  problem). 
Let  vp  be  an  approximation  to  the  solution  vector  and  define  the  discrete  residual,  rp(vp), 
by 


rp(vp)  =  fp  —  Rp(vp). 

In  a  basic  two-level  multigrid  method,  the  exact  solution  on  a  coarse  level  is  used  to  correct 
the  solution  on  the  fine  level.  This  multigrid  scheme  is  given  as  follows: 

•  Perform  v\  smoothing  alterations  on  the  fine  level:  vp,n+1  =  vp,n  —  P_1rp(vp,n) 
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•  Restrict  the  state  and  residual  to  the  coarse  level:  Vg_1  =  Ip_1vp,  rp_1  =  Ip  Vp. 

•  Solve  the  coarse  level  problem:  Rp_1(vp_1)  =  Rp_1(vg_1)  +  rp_1. 

•  Prolongate  the  coarse  level  error  and  correct  the  fine  level  state:  vp  =  vp+Ip_1(vp_1  — 

•  Perform  zv2  smoothing  interations  on  the  fine  level:  vp,n+1  =  vp,n  —  P_1rp(vp,n). 

Ip ~1  is  the  residual  restriction  operator,  and  /^_1  is  the  state  prolongation  operator.  Ip  1 
is  the  state  restriction  operator  and  is  not  necessarily  the  same  as  the  residual  restriction. 
Alternatively,  the  FAS  coarse  level  equation  can  be  written  as 

RP-I^P-1)  =  JP-lfP  +  7-P-l, 

r^1  =  R p-1(/p-1vp)-/p-1Rp(vp). 

The  first  equation  differs  from  the  original  coarse  level  equation  by  the  presence  of  the  term 
Tp_1,  which  improves  the  correction  property  of  the  coarse  level.  In  particular,  if  the  fine 
level  residual  is  zero,  the  coarse  level  correction  is  zero  since  vp_1  =  1 . 

3.3.3  V-cycles  and  FMG 

To  make  multigrid  practical,  the  two-level  correction  scheme  is  extended  to  V-cycles  and  Full 
Multigrid  (FMG).  In  a  V-cycle,  one  or  more  levels  are  used  to  correct  the  fine  level  solution. 
Descending  from  the  finest  level,  after  restriction,  u\  smoothing  steps  are  performed  at  each 
level  until  the  coarsest  level  is  reached.  On  the  coarsest  level,  the  problem  may  be  solved 
exactly  or  smoothed  a  relatively  large  number  of  times.  Ascending,  the  problem  is  smoothed 
z/2  times  on  each  level  after  prolongation  until  the  finest  level  is  reached.  This  procedure 
constitutes  one  V-cycle  which  is  also  refered  to  as  one  multigrid  iteration. 

Using  only  V-cycles  to  obtain  high-order  solutions  is  impractical  because  it  requires 
starting  the  calculation  on  the  finest  level,  where  there  are  the  most  degrees  of  freedom 
and  smoothing  is  most  expensive.  This  problem  can  be  eliminated  by  using  FMG,  where 
the  calculation  is  started  on  the  coarsest  level.  After  converging  or  partially  converging 
the  solution,  it  is  prolongated  to  the  next  finer  level.  Running  V-cycles  at  this  level,  the 
solution  is  partially  converged  and  then  prolongated  to  the  next  finer  level.  This  procedure 
continues  until  the  desired  solution  order  is  reached. 

For  more  details  on  the  prolongation  and  restriction  operators  or  the  multigrid  imple¬ 
mentation,  see  [22,  23]. 
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Chapter  4 


Stability  Analysis 


To  determine  the  stability  of  the  smoothers  discussed  in  Chapter  3,  Fourier  (Von  Neumann) 
analysis  is  performed  for  convection-diffusion  in  one  and  two  dimensions  with  periodic 
boundary  conditions.  The  analysis  follows  that  of  [22]  for  advection.  The  convection- 
diffusion  problem  is  given  by 

V-Vu-vV2u  =  /(f), 

aux  —  vuxx  =  f(x)  on  [-1,1], 
aux  +  buy  —  v(uxx  +  uyy)  =  f  (x,  y)  on  [-1, 1]  x  [-1, 1], 

For  this  analysis,  the  velocity,  V,  is  constant,  u  is  the  concentration  variable,  and  /  is  a 
source  term. 


4.1  Outline  of  Analysis 

To  begin,  the  domain  is  triangulated,  and  the  convection-diffusion  equation  is  discretized 
following  the  steps  in  Chapter  2.  The  resulting  discrete  system  is  linear  and  will  be  written 
Au  =  f,  where  u  is  the  exact  solution.  Denoting  the  current  solution  guess  as  vn,  the 
general  iterative  solution  procedure  as  defined  in  Chapter  3  is  given  by 

vn+1  =  vn  -  P-V, 

where 

rn  =  Av"  —  f .  (4.1) 
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Defining  the  error  at  iteration  n  to  be  e"  =  v”  —  u,  Eqn.  4.1  can  be  written  in  terms  of  the 
error: 

Aen  =  rn.  (4.2) 

Thus,  the  iterative  scheme  can  also  be  put  in  terms  of  the  error: 

en+i  =  Sen; 

where  S  =  I  —  P_1A  is  the  iteration  matrix  and  I  is  the  identity  matrix.  The  spectral 
radius  of  the  iteration  matrix,  |/o(S)|,  determines  the  growth  or  decay  of  the  error.  Thus, 
to  determine  the  stability  of  the  iterative  scheme,  it  is  necessary  to  compute  the  eigenvalue 
footprint  of  this  matrix.  For  stability,  the  eigenvalues  of  S  must  lie  in  the  unit  circle  centered 
at  the  origin.  This  stability  condition  requires  that  the  eigenvalues  of  — P_1A  lie  in  the 
unit  circle  centered  at  (—1,0).  Thus,  eigenvalue  footprints  of  —  P-1A  are  computed  for 
both  element  and  element-line  Jacobi  relaxation  via  Fourier  analysis.  The  specifics  of  the 
discretization  and  Fourier  analysis,  including  results,  are  presented  in  Sections  4.2  and  4.3. 


4.2  One  Dimensional  Analysis 


To  discretize,  the  domain,  [—1, 1],  is  divided  into  a  triangulation,  T)( ,  of  N  elements,  k.  of 
size  Ax  =  2/A,  such  that  (JKeT  n  =  [—1,1].  Then,  using  the  solution  and  test  function 
space  V?  defined  by  Eqn.  2.2,  the  DG  discretization  of  Eqn.  4.1  takes  the  following  form: 
find  Uh  £  such  that  Mvh  G  V/, 


Y,  'H(uh)v£\KR-Tt{uh)v£ \KL-  /  auhvh)Xdx+  /  vuhjXvKxdx 

„  rji  _  J  K,  K, 


*eTh 


Y  [HI  +  lUhi{"Vh,x}\  =  /  vhfdxi  (4-3) 

r  K€Th 


where  kl  and  hr  denote  the  left  and  right  boundaries,  respectively,  of  the  element  k.  and 
r  is  the  union  of  element  boundaries  over  the  entire  domain.  Full-upwinding  is  used  for  the 
inter-element  inviscid  flux,  such  that  TL{uh)  is  given  by, 

'H{uh)  =  \{a-  | a\)uh,R  +  ^(a  +  | a\)uh,L, 


where  Uh,L  and  Uh,R  are  values  of  Uh  taken  from  the  left  and  right  elements  at  an  interface. 
In  1-D,  the  local  lifting  operator  is  a  scalar  defined  by  the  following:  find  i5f  £  such  that 
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Vrh  G  V?, 


[  Th5LJR  dx  =  -  -  [vTh{uh,L  -  uhiR)\a  , 

Jkl/R  Z  1 


where  ct/  denotes  a  single  inter-element  boundary,  and  nL/R  represents  the  elements  to  the 
left  and  right  of  this  boundary.  Finally,  for  uniform  spacing,  r]f  =  2  at  each  interface. 


Using  a  basis  {4>j}  for  the  finite  element  space  Vf:,  the  concentration,  Uh ,  is  represented 
by  Uh(x)  =  Y2j  Uj<frj(x).  Thus,  Eqn.  4.3  can  be  written  concisely  as  Au  =  f,  where  u  is  the 
vector  of  solution  coefficients,  Uj.  For  this  analysis,  standard  Lagrange  basis  functions  are 
used. 


Assuming  the  error  varies  sinusoidally  on  the  elements,  it  can  be  decomposed  into  N 
modes, 

N/2 

e”  =  £  e"(03), 

j=-N/  2+1 

where  the  jth  error  mode  is  given  by, 


e”(%)  = 


e”(% 


L  J 


e'r(dj)  =  vn(^)exp  (irdj), 


(4.4) 


and  Oj  =  j  tv  Ax  and  r  is  the  element  index.  Thus,  e™(0j)  is  a  vector  of  length  p  +  1 
corresponding  to  the  jth  error  mode  on  element  r. 


Using  the  fact  that  the  stencil  is  element  compact,  the  Eqn.  4.2  can  be  written,  for  any 
element  r,  as 


:w^n 


U  +  A°e?  +  A 


'  E-n 


■'r+l 


=  r„ 


(4.5) 


where  Aw ,  A0,  and  AE  are  the  (r,r  —  1),  (r,  r),  and  (r,  r  +  1)  blocks  of  the  matrix  A. 
Because  the  boundary  conditions  are  periodic,  if  r  =  1,  r  —  1  refers  to  element  N,  and  if 
r  =  N,  r  +  1  refers  to  element  1.  Substituting  the  form  of  the  error  from  Eqn.  4.4  into 
Eqn.  4.5  gives 


A11  exp(— iOj)  +  A0  +  Ae  exp (i6j) 


en 


=  r 


n 
r  • 


Thus,  the  system  of  N(p  +  1)  equations  governing  the  error  represented  by  Eqn.  4.2  can  be 
reduced  to  a  system  of  (p  +  1)  equations  for  each  error  mode.  Furthermore,  the  iterative 
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scheme  can  also  be  reduced  such  that  the  relaxation  of  the  jth  error  mode  is  given  by 

e^+1(9j)  =  S (6j)e^(9j)  =  Sn+1  (9j)v° (9j)  exp(ir9j) , 

where  S (9j)  is  a  (p  +  1)  x  (p  +  1)  matrix  corresponding  to  the  iteration  matrix,  S,  for 
sinusoidal  error  variation.  Thus,  to  determine  the  stability  of  an  iterative  scheme,  one  must 
compute  the  eigenvalues  of  S (9j)  for  all  j. 

For  the  element  Jacobi  smoother,  the  preconditioner,  P,  is  the  block  diagonal  of  A. 
Thus,  for  sinusoidal  error,  the  (p+  1 )  x  (p+  1)  equivalent  of  the  iteration  matrix  is  given  by 

S(9J)  =  I-P~1A(0j), 


where 

P-1A(^-)  =  (A°)-1(Awexp(— +  A0  +  AE  exp(i9j)). 

Footprints  of  —  P_1A  for  the  1-D  element  Jacobi  smoother  are  shown  in  Figure  4-1.  In 
1-D,  the  footprints  depend  only  on  the  solution  order,  p,  and  the  element  Reynolds  number, 
Re  =  aAx/v.  The  figures  show  footprints  for  p  =  0, 1,2,3  at  four  Re.  Note  that  all  the 
eigenvalues  are  stable  and  that,  as  Re  — >  oo,  all  the  eigenvalues  associated  with  p  >  0  are 
centered  at  the  origin  [22]. 

In  1-D  the  element-line  Jacobi  smoother  becomes  an  exact  solve.  Thus,  this  smoother 
is  only  examined  in  2-D. 


4.3  Two  Dimensional  Analysis 

In  2-D,  the  domain  is  subdivided  into  NxNy  rectangular  elements,  k.  of  size  Ax  by  Ay  where 
Ax  =  2/Nx  and  Ay  =  2 /Ny.  The  discretization  procedure  is  analogous  to  that  for  the  1-D 
case.  The  2-D  basis  is  given  by  the  tensor  product  of  the  1-D  basis:  4>ag(x,  y)  =  (j)Q(x)(j)p(y), 
where  (j>a  and  4>g  are  1-D  Lagrange  basis  functions.  Thus,  there  are  (p+  l)2  basis  functions 
per  element. 

Indexing  elements  by  the  ordered  pair  (r,  s),  error  modes  have  the  form 


er,s(^j^k)  =  vn  exp  (ir9j  +  is9k ), 


where  9j  =  jirAx,  9k  =  kirAy,  and  j,  k  £  (—IV/ 2  +  1, ...,  A/2).  Thus,  for  the  element  (r,  s ), 
Eqn.  4.2  becomes 


An/  exp(— i9j)  +  As  exp {—i9k)  +  A0  +  AE  exp (i9j)  +  AiV  exp (i9k) 


=  r 


n 

r,s’ 
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(c)  Re  =  1000 


(d)  Re  =  10000 


Figure  4-1:  Eigenvalue  footprints  for  element  Jacobi  preconditioned  1-D  convection-diffusion 


Figure  4-2:  Stencil  of  element  0 


where  each  A  matrix  is  size  (p+  l)2  x  (p+  l)2. 

Thus,  the  (p  +  l)2  X  (p  +  l)2  iteration  matrix  corresponding  to  S  for  sinusoidal  error  is 
given  by  the  following: 

s(0j,ek)  =  i-P-1(ej,ek)A(ej,ek). 


The  preconditioner  in  the  element  smoothing  case  is  still  simply  the  block  diagonal: 


P  =  A0. 


The  form  of  the  element-line  smoother  depends  on  the  direction  of  the  lines.  As  de¬ 
scribed  in  Section  3.2,  lines  are  formed  based  on  the  coupling  between  elements  for  a  p  =  0 
discretization  of  convection-diffusion.  Figure  4-2  shows  the  dependence  of  ro  on  the  sur¬ 
rounding  elements.  For  p  =  0,  the  residual  and  error  on  each  element  are  constant,  thus, 
the  residual  on  element  zero  is 


e_B 


ro  = 


<9ro  dro  <9ro  <9ro  <9ro 
deE  dew  deo  dew  des 


eN 

eo  j 
ew 


es 


where 


dip 

(^Ay(a-\a\)  - 

/  A  \ 

des 

dip 

deN 

(jAx(6-|6|) 

dip 
de  o 

= 

(As/|o|  +  +  2i/(£§  +  If)) 

dip 
de\ y 

l-lAi/(o  +  |o|)  - 

dip 

des 

(-iAx(6+|6|)-|f„) 
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Note  that  because  the  problem  has  constant  coefficients  and  the  grid  is  uniform, 


drE 

dr0 

de0 

dew 

drN 

dr0 

de0 

des ’ 

dr  w 

dr0 

de0 

deE’ 

9rs 

dr0 

de0 

deN 

These  relationships  imply  that  the  vertical  face  and  horizontal  face  connectivity  values  are 
constant  throughout  the  domain.  Thus,  if  the  horizontal  connectivity  is  greater  than  the 
vertical  connectivity, 


max 


( 

9r0 

dr0 

l  A>  m  o  v  1 

<9r0 

dro 

V 

deE 

1 

dew 

1  /  iildA.  1 

deN 

5 

des 

which  simplifies  to 

Ay  (|o1  +  £)  2  (|61  +  iy)  ’ 

the  lines  are  horizontal  everywhere.  Otherwise,  they  are  vertical. 
For  the  horizontal  case,  the  preconditioner  is 


P (6j)  =  A1'  exp(— idj)  +  A0  +  AE  exp(i6j), 


and  for  the  vertical  case, 

P {Ok)  =  As  exp (—iOk)  +  A0  +  AiV  exp iiOk). 

In  2-D,  the  eigenvalue  footprints  depend  of  five  parameters:  smoother  choice;  solution 
order,  p:  element  Reynolds  number,  Re  =  aAx/w,  element  aspect  ratio,  AR  =  Ax/ A y\  and 
flow  angle,  tana  =  b/a.  Footprints  of  —  P-1A  for  both  element  and  element-line  Jacobi 
smoothing  are  shown  in  Figure  4-3  and  Figure  4-4.  Two  flow  conditions  are  shown:  a  low 
AR,  low  Re  case,  and  a  high  AR,  high  Re  case.  The  cases  are  meant  to  be  representative 
of  the  range  of  flow  conditions  one  would  encounter  in  a  typical  aerodynamic  simulation. 
The  low  AR,  low  Re  case  approximates  the  conditions  one  would  expect  near  the  leading 
edge  of  an  airfoil,  while  the  high  AR,  high  Re  case  is  representative  of  a  boundary  layer  in  a 
high  speed  flow.  While  only  two  cases  are  shown,  a  comprehensive  study  of  the  parameter 
space  was  conducted  using  1  <  AR  <  500,  0  <  Re  <  10000,  and  1°  <  a  <  25°.  As  in  the 
1-D  case,  the  smoothers  are  stable  independent  of  p  and  flow  condition. 
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Figure  4-3:  Eigenvalue  footprints  for  element  Jacobi  preconditioned  2-D  convection-diffusion 


p  =  o 


p  =  2 


(a)  Horizontal  lines,  AR  =  1,  Re  =  10,  a  =  1°  (b)  Vertical  lines,  AR  =  100,  Re  =  10000,  a  =  1° 


Figure  4-4:  Eigenvalue  footprints  for  element-line  Jacobi  preconditioned  2-D  convection- 
diffusion 


Chapter  5 


Numerical  Results 


In  this  chapter,  three  test  cases  are  considered.  For  the  first  two  cases,  source  terms  are 
added  to  the  compressible  Navier-Stokes  equations  such  that  a  desired  analytic  function  is 
the  exact  solution.  Thus,  the  L2-norm  of  the  error  in  the  solution  may  be  computed.  Mesh 
refinement  studies  for  these  cases  show  optimal  order  of  accuracy. 

The  last  case  presented  is  that  of  M  =  0.5,  Re  =  5000,  laminar  flow  over  a  NACA 
0012  airfoil  at  a  =  0°.  A  grid  refinement  study  shows  that  the  drag  on  the  airfoil  can  be 
computed  accurately  in  less  time  using  high-order  rather  than  highly  refined  meshes. 


5.1  Poiseuille  Flow 

To  numerically  verify  the  order  of  accuracy  of  the  discretization,  the  first  test  case  considered 
is  that  of  fully  developed  flow  in  a  straight  channel.  Thus,  the  exact  solution  is  of  the  form: 

p  =  constant, 

v  =  0. 

Of  course,  the  classical,  Poiseuille  form  is  a  solution  of  the  incompressible  Navier-Stokes 
equations.  To  make  this  form  satisfy  the  compressible  Navier-Stokes  equations,  constant 
viscosity  is  assumed,  and  a  source  term  is  added  to  the  right-hand  side  of  the  energy 
equation.  Thus,  the  equations  are 


u4  +  V  •  ^i(u)  -  V  •  ^(u,  Vu)  =  S, 


(5.1) 
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where, 


S  = 


0 

0 

0 


V  b  (j&)  y^y  ~ 6)  [^=r  -  h(y  - 6)] 

Finally,  in  terms  of  the  conservative  state  variables,  the  exact  solution  is  given  by 


/  0  \ 


u  = 


P 

pu 

pv 


(  1  \ 

b  (I?)  y(y~b ) 


\PE  )  \  (Po  +  £x)  +  U2  J 


For  this  test,  the  parameters  p,  po,  and  ^ 


are  chosen  as 


p  =  lx  1(T4, 


Po 

dp 

dx 


1, 

-8  x  1(T4. 


The  domain  of  interest  is  a  1  x  1  square  shown  in  Figure  5- 1(a).  On  the  left  and 
right  boundaries,  the  boundary  state  is  set  to  the  exact  solution.  On  the  top  and  bottom 
boundaries,  no  slip,  isothermal  wall  conditions  are  imposed. 


(a)  Dimensions  and  BCs 


(b)  Medium  mesh,  176  elements 


Figure  5-1:  Domain  for  Poiseuille  flow  test  case 
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(a)  Error  vs.  Grid  Spacing 


(b)  Error  vs.  DOF 


Figure  5-2:  Accuracy  results  for  0  <  p  <  3  for  Poiseuille  flow  test  case 

The  coarsest  grid  used  contains  44  elements.  The  next  finer  grid  (176  elements),  which 
is  shown  in  Figure  5-l(b),  was  generated  by  subdividing  each  element  into  four  elements 
such  that  the  grid  spacing,  h,  is  halved.  This  procedure  was  repeated  a  total  of  three  times, 
resulting  in  four  meshes  containing  44,  176,  704,  and  2816  elements.  Note  that,  since  the 
channel  has  straight  geometry,  it  can  be  represented  exactly  with  straight-edged  ( q  =  1) 
triangles.  From  this  point  forward,  q  is  used  to  represent  the  order  of  the  basis  used  for  the 
geometry  interpolation. 

Results  of  the  mesh  refinement  study  are  shown  in  Figure  5-2.  Figure  5-2(a)  shows 
that  the  L2-norm  of  the  error  is  converging  at  approximately  0{hp+1).  Furthermore,  as 
demostrated  by  Figure  5-2(b),  the  higher-order  solutions  are  dramatically  more  accurate 
per  degree  of  freedom  (DOF).  For  the  finest  mesh,  p  =  1  solution,  which  has  3.38  x  104 
DOF,  the  error  is  approximately  4.27  x  10  ,  while  the  coarsest  mesh,  p  =  3  solution  has 

only  7.04  x  103  DOF  and  gives  an  error  of  2.36  x  10-4. 


5.2  Circular  Channel  Flow 

The  goal  of  the  second  test  case  was  to  numerically  verify  that  the  discretization  achieves 
optimal  order  of  accuracy  for  a  problem  with  curved  boundaries.  The  domain,  shown  in 
Figure  5-3(a),  is  a  180°  section  of  a  circular  channel.  Figure  5-3(b)  shows  the  coarsest  grid 
used,  which  contains  100  elements.  Note  that  on  the  circular  boundaries  the  geometry  is 
represented  using  a  q  =  2  Lagrange  basis  where  the  midpoint  node  on  each  boundary  edge 
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r  o  =  2.0 

(a)  Dimensions  and  boundary  conditions 


(b)  Coarse  mesh,  100  elements,  q  =  2 


Figure  5-3:  Domain  for  circular  channel  test  case 
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(a)  Error  vs.  Grid  Spacing  (b)  Error  vs.  DOF 


Figure  5-4:  Accuracy  results  for  0  <  p  <  3  for  circular  channel  test  case 
is  located  on  the  actual  geometry. 

As  in  Section  5.1,  a  source  term  is  added  to  the  right-hand  side  of  the  compressible 
Navier-Stokes  equations  such  that  the  exact  solution  to  the  equations  is 


p  \ 

(  1 

pu 

(r  —  ri){r  —  r0)  sin  9 

pv 

(r  —  ri)(r  —  r0)  cos  8 

V  ()E  ) 

\  5+  l(r  -  n)2(r  -  ra)2  y 

Figures  5-4(a)  and  5-4(b)  show  the  Z^morm  of  the  error  versus  number  of  elements  and 
DOF.  As  for  the  the  straight  channel  Poiseuille  flow  problem,  optimal  order  of  accuracy  is 
observed,  and  the  higher-order  solutions  are  seen  to  provide  significantly  lower  error  with 
fewer  DOF.  For  p  =  3,  only  4000  DOF  are  required  to  achieve  an  error  of  1.37  x  10-4,  while 
for  p  =  1,  the  error  is  2.64  x  10-4  for  7.68  x  104  DOF. 

While  conducting  this  test,  an  unexpected  phenomenon  was  observed.  When  using  q  =  3 
geometry  representation,  the  order  of  accuracy  dropped  from  0(hp+1 )  to  0{hP+ 1//2)  for  p  =  2 
and  p  =  3.  While  it  is  not  clear  at  this  time,  this  problem  may  be  the  result  of  oscillations  in 
the  geometry  interpolation.  Figure  5-5  shows  the  radius  error — the  difference  between  the 
interpolated  radius  and  the  exact  radius — for  a  coarse  grid  element  on  the  inner  wall.  The 
radius  error  is  plotted  versus  the  location,  er,  on  the  reference  edge.  Clearly,  both  the  q  =  2 
and  q  =  3  interpolations  produce  oscillatory  geometry.  In  fact,  the  q  =  2  oscillations  are 
larger  in  amplitude.  Yet,  with  q  =  2  iterpolants,  optimal  accuracy  is  achieved.  Further  work 
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Figure  5-5:  Radius  error  for  inner  wall,  coarse  grid  boundary  face 

in  this  area  is  necessary  to  fully  understand  the  requirements  on  the  geometry  convergence 
to  achieve  optimal  order  of  accuracy.  Specifically,  it  is  unclear  not  only  what  norm  of  the 
geometry  error  is  appropriate  but  also  what  rate  this  norm  should  converge  to  guarantee 
0(hp+1)  convergence  of  the  solution  error. 


5.3  NACA  0012  Airfoil 

The  last  test  case  is  that  of  M  =  0.5,  Re  =  5000,  laminar  flow  over  a  NACA  0012  airfoil  at 
a  =  0°.  Three  sets  of  results  are  shown:  accuracy  results,  iterative  convergence  results,  and 
timing  results.  A  set  of  four  structured  grids  was  used  to  generate  the  data.  The  meshes 
were  created  by  modifying  a  baseline  grid  provided  by  Swanson  [43] . 

The  leading  edge  of  the  original  Swanson  mesh  was  refined  by  inserting  points  along 
the  airfoil.  The  points  were  used  to  create  new  elements  having  lower  aspect  ratios  than 
the  original  elements  such  that,  when  the  boundary  elements  were  converted  from  q  =  1 
to  q  =  3  by  moving  the  higher-order  points  to  the  true  airfoil  surface,  no  negative  area 
elements  were  created.  Then,  to  create  a  family  of  meshes  with  varying  h,  the  mesh  was 
coarsened  by  removing  every  other  node  in  both  the  streamwise  and  normal  directions.  The 
coarsening  procedure  was  repeated  three  times,  resulting  in  four,  nested  meshes  containing 
672,  2688,  10752,  and  43008  elements.  The  coarsest  mesh  is  shown  in  Figure  5-6.  On 
the  airfoil  boundary,  q  =  3  elements  were  used.  The  boundary  nodes  were  placed  using  a 
modified  version  of  the  analytic  definition  of  the  NACA  0012: 

y  =  ±0.6(0.2969Vx  -  0.1260x  -  0.3516x2  +  0.2843x3  -  0.1036x4).  (5.2) 
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Figure  5-6:  Coarse  NACA  0012  grid,  672  elements 

In  Eqn.  5.2,  the  x4  coefficient  of  the  traditional  definition  of  the  NACA  0012  has  been 
modified  such  that  the  trailing  edge  (x  =  1)  has  zero  thickness. 

On  the  airfoil,  the  no  slip,  adiabatic  boundary  condition  is  used.  At  the  exit  plane, 
the  subsonic  outflow  condition  is  imposed,  thus,  the  pressure  is  set  to  the  freestream  value. 
Along  the  rest  of  the  boundary,  the  full  state  condition  is  used,  where  the  boundary  state 
is  set  to  the  freestream  state. 

Unless  otherwise  noted,  all  results  were  obtained  using  the  multigrid  solver  with  element¬ 
line  Jacobi  smoothing.  All  computations  were  initialized  from  a  fully  converged  p  =  0 
solution.  When  starting  higher-order  calculations,  five  V-cycles  on  each  level  were  run 
before  prolongating  to  the  next  finer  level,  and  each  V-cycle  contained  four  pre-  and  post¬ 
smoothing  iterations  with  100  smoothing  iterations  on  the  coarsest  level  (p  =  0). 

5.3.1  Accuracy  Results 

The  computed  drag  for  each  grid  and  interpolation  order  is  shown  in  Figure  5-7.  Figure  5-8 
shows  the  drag  versus  DOF.  For  comparison,  results  computed  by  FUN2D,  which  uses  a 
node-centered  finite  volume  algorithm,  on  the  same  meshes  are  shown.  While  it  appears 
that  FUN2D  may  be  converging  to  a  different  final  answer  for  the  drag,  qualitatively,  the 
convergence  of  the  FUN2D  and  p  =  1  methods  in  terms  of  DOF  is  very  similar. 

However,  the  p  =  2  and  p  =  3  results  are  dramatically  better.  On  the  2688  element 
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(a)  All  results 


(b)  Close-up  of  (a) 


Figure  5-7:  Drag  versus  number  of  elements 
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Figure  5-8:  Drag  versus  degrees  of  freedom 
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grid,  both  the  p  =  2  (6.45  x  104  DOF)  and  p  =  3  (1.08  x  105  DOF)  drag  results  are  within 
0.1  drag  counts  of  the  finest  grid,  p  =  3  result.  Alternatively,  the  drag  error  for  p  =  1  on 
the  43008  element  mesh  (5.16  x  105  DOF)  is  0.3  drag  counts.  Thus,  high-order  solutions 
provide  significantly  better  accuracy  in  fewer  DOF  than  p  =  1  solutions  on  very  refined 
meshes. 

It  should  be  noted  that  in  the  results  presented  here,  the  drag  error  is  not  converging  at 
0(hp+1).  While  the  p  =  1  results  are  close — the  convergence  rate  based  on  the  p  =  1,  10752 
and  43008  element  grids  is  1.78 — the  p  =  2  and  p  =  3  results  are  not.  The  loss  of  optimal 
order  of  accuracy  could  be  caused  by  many  factors.  First,  as  noted  in  Section  5.2,  the 
discretization  is  very  sensitive  to  geometry  errors.  While  a  q  =  3  geometry  representation 
was  used,  the  interpolated  geometry  is  oscillatory,  and  it  is  possible  that  the  geometry 
errors  converging  at  suboptimal  order  of  accuracy  could  be  affecting  the  order  of  accuracy 
of  the  solution.  Second,  the  smoothness  of  the  exact  solution  could  be  affecting  the  order 
of  accuracy.  As  is  well  known,  the  exact  solution  has  a  singularity  at  the  trailing  edge.  No 
effort  has  been  made  here  to  over-refine  the  trailing  edge  region,  and,  thus,  this  singularity 
may  be  adversely  affecting  the  convergence  rate. 

While  the  above  factors  appear  to  be  the  most  likely  culprits,  many  other  problems 
could  be  degrading  the  order  of  accuracy.  For  example,  the  boundary  conditions  have  not 
been  thoroughly  investigated,  and  the  grid  spacing  does  not  vary  smoothly  throughout  the 
domain.  Finally,  it  should  be  noted  that,  while  optimal  order  of  accuracy  has  been  proven 
for  DG  for  linear  hyperbolic  and  linear  elliptic  problems,  no  proof  currently  exists  for  the 
Navier-Stokes  equations. 

5.3.2  Iterative  Convergence  Results 

Iterative  convergence  results  are  shown  in  Figures  5-9  and  5-10.  Residual  convergence  for 
1  <  p  <  3  are  shown  for  three  meshes  in  Figure  5-9.  The  spikes  clearly  observable  for  p  =  2 
at  five  iterations  and  for  p  =  3  at  both  five  and  ten  iterations  are  where  the  solver  has 
begun  to  run  V-cycles  on  the  next  finer  level. 

As  demonstrated  in  [22]  for  the  Euler  equations,  the  asymptotic  convergence  rate  of  the 
residual  is  approximately  independent  of  p.  However,  a  slight  h  dependence  is  observed. 
Thus,  while  the  residual  on  the  672  element  mesh  is  fully  converged  in  less  than  50  iterations, 
the  residual  is  approximately  10-8  after  50  iterations  on  the  10752  element  mesh. 

Figure  5-10  shows  the  drag  history.  At  each  iteration,  the  iterative  drag  error — the 
difference  between  the  drag  at  that  iteration  and  the  final  converged  drag  for  that  case — is 
plotted.  Note  that,  unlike  the  residual  convergence  histories,  only  V-cycles  on  the  finest 
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(a)  672  Elements 


(b)  2688  Elements 


(c)  10752  Elements 


Figure  5-9:  Residual  convergence  history 
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level  are  shown. 


The  results  show  that  the  drag  converges  to  engineering  tolerance  much  faster  than  the 
residual  converges  to  zero.  In  fact,  the  iterative  drag  error  usually  converges  to  less  than 
0.01  drag  counts  in  less  than  10  multigrid  iterations  on  the  finest  level,  and  more  than  15 
are  never  required.  In  practice,  this  means  that  the  calculation  may  be  stopped  before  the 
residual  has  converged  to  zero. 


(a)  672  Elements 


(b)  2688  Elements 


(c)  10752  Elements 


Figure  5-10:  Drag  convergence  history 
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CPU  Time  (s)  CPU  Time  (s) 


(a)  All  results 


(b)  Close-up  of  (a) 


Figure  5-11:  Absolute  drag  versus  CPU  time 

5.3.3  Timing  Results 

Two  sets  of  timing  results  are  presented:  a  comparison  of  the  time  required  to  compute 
the  drag  to  engineering  tolerance  using  multigrid  with  element-line  Jacobi  smoothing  for 
varying  p  and  h  and  a  comparison  of  the  time  required  to  drive  the  residual  to  zero  using 
the  element  Jacobi  smoother,  element-line  Jacobi  smoother,  multigrid  with  element  Jacobi 
smoothing,  and  multigrid  with  element-line  Jacobi  smoothing.  All  timing  study  calculations 
were  run  on  a  single  3.0  gigahertz  Intel  Pentium  4  processor  with  896  megabytes  of  memory. 

Figure  5-11  shows  the  time  required  to  compute  the  drag  using  multigrid  with  element¬ 
line  Jacobi  smoothing  for  varying  p  and  h,  and  Figure  5-12  shows  the  drag  error  versus 
time,  where  drag  error  is  the  difference  between  the  drag  for  a  given  case  and  the  p  = 
3  result  on  the  43008  element  grid.  The  values  plotted  were  obtained  using  a  drag 
based  stopping  criterion.  Namely,  when  the  change  in  the  drag  for  one  iteration  was  less 
that  0.01  drag  counts,  10000|c^+1  —  c^|  <  0.01,  for  four  iterations,  the  calculation  was 
stopped  and  the  drag  and  time  were  plotted.  This  stopping  criterion  is  intended  to  minimize 
the  computational  time  required  to  compute  the  drag  given  the  knowledge  that  the  drag 
converges  to  engineering  tolerance  faster  than  the  residual  converges  to  zero.  Furthermore, 
it  is  intended  to  simulate  engineering  practices  in  that  it  does  not  require  knowledge  of  the 
final,  converged  drag. 

Clearly,  the  higher-order  solutions  give  smaller  error  in  less  time  than  the  p  =  1  solution. 
It  takes  approximately  2.98  x  105  seconds  to  obtain  the  most  accurate  p  =  1  solution,  which 
has  an  error  of  0.32  drag  counts.  Alternatively,  the  coarsest  p  =  2  solution  has  an  error  of 
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Figure  5-12:  Drag  error  versus  CPU  time 

0.29  counts  and  was  obtained  in  only  548  seconds. 

Figure  5-13  compares  the  four  solvers  developed  in  this  work:  element  Jacobi,  element¬ 
line  Jacobi,  multigrid  with  element  Jacobi  smoothing,  and  multigrid  with  element-liue  Ja¬ 
cobi  smoothing.  As  expected,  the  figure  shows  that  the  incorporation  of  the  element-line 
Jacobi  smoother  significantly  improves  the  performance  of  the  solver.  Thus,  the  element-line 
Jacobi  solver  is  more  efficient  than  the  element  Jacobi  solver,  and  multigrid  with  element¬ 
line  Jacobi  smoothing  is  more  efficient  that  multigrid  with  element  Jacobi  smoothing  for 
all  p.  Furthermore,  in  all  cases,  multigrid  with  element-line  Jacobi  smoothing  is  the  best 
solver,  and,  as  p  increases,  it  becomes  relatively  more  efficient.  For  p  =  1,  the  multigrid 
with  element-line  Jacobi  smoothing  is  only  slightly  more  efficient  than  pure  element-line  Ja¬ 
cobi.  After  3000  seconds,  the  residual  is  approximately  one  order  of  magnitude  lower  using 
multigrid.  However,  for  p  =  3  the  element-line  Jacobi  and  multigrid  with  element  Jacobi 
schemes  both  reduce  the  residual  to  approximately  10~3  in  8000  seconds  while  multigrid 
with  element-line  Jacobi  drives  the  residual  to  10-9  in  the  same  time. 
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(a)  V  =  1 


(b)  V  =  2 


(c)  V  =  3 


Figure  5-13:  Solver  comparison  for  2688  element  mesh 
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Chapter  6 


Conclusions 


In  an  effort  to  advance  the  state  of  the  art  in  CFD  in  applied  aerodynamics,  a  high-order 
algorithm  for  the  compressible  Navier-Stokes  equations  was  presented.  The  algorithm  em¬ 
ploys  a  DG  discretization  with  an  element-compact  stencil  and  a  p-multigrid  solver  with 
element-line  Jacobi  smoothing. 

The  discretization  has  been  shown  to  achieve  optimal  order  of  accuracy  for  simple  test 
problems.  While  optimal  order  of  accuracy  was  not  obtained  for  the  more  practical  airfoil 
problem,  the  higher-order  discretizations  were  shown  to  provide  significantly  less  error  with 
fewer  elements  and  DOF  than  second-order  discretizations.  Furthermore,  timing  studies 
showed  that  the  reductions  in  elements  and  DOF  translated  into  savings  in  computational 
time  required  to  compute  the  drag  to  within  0.5  drag  counts  of  an  order  of  magnitude  or 
more. 

Via  Fourier  analysis  of  convection-diffusion,  the  element  and  element-line  Jacobi  smoothers 
were  shown  to  be  stable  independent  of  order.  Furthermore,  the  residual  convergence  rate 
was  observed  to  be  independent  of  order  and  only  weakly  dependent  on  grid  spacing  when 
using  multigrid  with  element-line  Jacobi  smoothing.  A  timing  study  on  the  NACA  0012 
test  case  using  all  four  solvers  confirmed  that,  especially  for  higher-order  discretizations, 
multigrid  with  element-line  Jacobi  smoothing  is  the  most  efficient  solution  technique. 

While  these  results  are  encouraging,  much  work  remains.  First,  the  2-D  viscous  dis¬ 
cretization  must  be  extended  to  the  3-D  case.  Second,  the  geometry  representation  issue 
highlighted  by  the  circular  cylinder  test  case  must  be  resolved.  To  obtain  high-order  ac¬ 
curacy  for  general  problems,  a  complete  understanding  of  the  geometry  representation  and 
gridding  requirements  is  crucial.  Third,  turbulence  modeling  and  limiting  capabilities  must 
be  added  to  allow  consideration  of  a  wider  range  of  flows  of  engineering  interest.  In  terms 
of  turbulence  modeling,  the  first  model  considered  will  be  Spalart-Allmaras  as  this  model 
is  widely  used  in  aerospace  applications.  In  terms  of  limiting,  although  it  is  not  yet  clear, 
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the  most  attractive  approach  may  be  /rp-adaptivity  where  in  the  vicinity  of  a  shock  h  and 
p  are  lowered  to  avoid  oscillations  in  the  solution.  One  hurdle  in  the  introduction  of  such  a 
limiter  is  the  development  of  a  robust  shock  detection  scheme  such  that  h  and  p  are  lowered 
near  the  shock  but  high-order  accuracy  is  maintained  in  smooth  regions. 

Finally,  optimization  is  required  to  make  the  method  practical.  This  work  will  entail 
not  only  optimization  of  the  current  implementation  but  possibly  also  the  introduction  of 
approximations  to  lessen  the  impact  of  expensive  computations.  For  example,  in  the  current 
implementation,  the  block  diagonal  of  the  Jacobian  is  inverted  for  every  element  at  every 
iteration.  Especially  for  higher-order,  3-D  calculations,  this  inversion  requires  significant 
computational  work.  Thus,  if  an  approximate  inverse  could  be  introduced  without  severly 
degrading  the  convergence  rate,  significant  savings  could  be  obtained.  Work  in  this  area 
has  already  begun. 
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